TI/Optymalizacja

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue May 21 11:34:40 2019

@author: tgub

Wylosuj 10000 liczb z rozkładu normalnego z parametrami loc=1.23 i scale=2.0. Do wylosowanych danych dopasuj rozkład normlany trzema metodami

    METODA 1 - stwórz histogram otrzymanych wartości, znormalizuj go i metodą najmniejszych kwadratów dopasauj gęstość rozkładu do histogramu
    METODA 2 - dopasuj gęstość rozkładu do wylosowanych danych metodą największej wiarygodności
    METODA 3 - z wylosowanych danych stwórz dystrybuantę empiryczną. Metodą najmniejszych kwadratów dopasuj dystrybuantę rozkładu normalnego do dystrybuanty empirycznej.

"""



import scipy.optimize as so
import scipy.stats as ss
import pylab as py
import numpy as np
N=10000
x=2*np.random.randn(N)-1.23

#suma kwadratów
def squares(a,func,xlist,ylist):
    return sum([(func(xlist[i],a)-ylist[i])**2 for i in range(len(xlist))])
#definiujemy -funkcje wiarygodnosci
def logL(a,f,xwekt):
    return -sum([np.log(f(x,a)) for x in xwekt])
#gestosc rozkładu normalnego
def rho_normalny(x,a): 
    return 1/((2*np.pi)**0.5)/abs(a[1])*np.exp(-((x-a[0])**2)/(2*(a[1])**2))

def F_normalny(x,a):
    return ss.norm.cdf(x,loc=a[0],scale=a[1])

#METODA 1 - Dopasowanie metoda najmniejszych kwadratow do histogramu
 
#tworzymy histogram
hist,bins= np.histogram(x,bins=61)
#dlugosc przedzialu histogramowania
przedzial=bins[1]-bins[0]
#normalizujemy histogram aby moc go porownac z gestoscia
hist=hist*1.0/len(x)/przedzial
#liczymy wsp. srodkow przedzialow histogramowania
xhist=bins[:-1]+0.5*przedzial
#szukamy minimum funkcja fmin
fit1=so.minimize(squares,(0,1),args=(rho_normalny,xhist,hist))
print('wynik metody1 to '+str(fit1.x))
#ogladamy wynik
#py.plot(xhist,hist)
#xtest=np.linspace(-20,20,1001)
#ytest=[rho_normalny(a,*fit1.x) for a in xtest]
#py.plot(xtest,ytest)
#py.show()

#METODA 2 - Metoda najwiekszej wiarygodnosci

#szukamy minimum
fit2=so.minimize(logL,(0,1),args=(rho_normalny,x))
print('wynik metody2 to '+str(fit2.x))

#METODA 3 - dopasowanie dystrybuant
 
xx=sorted(x)
yy=np.linspace(0,1,N)
#szukamy minimum
fit3=so.minimize(squares,(0,1),args=(F_normalny,xx,yy))
print('wynik metody3 to '+str(fit3.x))
#ogladamy wynik
#cut=100
#py.plot(xx[cut:-cut],yy[cut:-cut])
#xtest=np.linspace(-20,20,1001)
#ytest=[F_normalny(x,fit3.x[0]) for x in xtest]
#py.plot(xtest,ytest)
#py.show()