FYSA1110 projektitehtävä

Sami Yrjänheikki

Tehtävänanto

Tässä projektissa simuloidaan normaalijakautuneiden muuttujien käyttäytymistä Pythonilla. Käytetään funktiota f(x, y) = 4xy + 3 ja selvitetään, kuinka normaalijakautuneiden muuttujien x ja y funktion arvot ovat jakautuneet.

Funktion jakauman teoreettinen ennustus

Tulosjakauman keskikohta saadaan suoraan muuttujien keskikohtien avulla siten, että μ = f(1.7, 2.6) = 4 1.7 2.6 + 3 = 20.68. Leveysparametri σ saadaan virheen etenemislaista seuraavasti:

yvel1.png

Kun käytetään muuttujien x ja y tilalla niiden odotusarvoja, saadaan

yvel2.png

Leveysparametri saa siis arvon σ ≈ 2.3.

Käytetyt parametrit

Määritellään ohjelman vaatimat parametrit. Tilastolliset parametrit on suoraan tehtävänannosta, loput ohjelman sisäisiä parametreja (ohjelma toimii hyvin oletusarvoilla - niitä ei tarvitse/kannata muuttaa).

In [19]:
# Tehtävässä annetut parametrit
mean_x = 1.7
sd_x = 0.2
mean_y = 2.6
sd_y = 0.15

# Tehtävässä annettu funktio f(x,y) = 4xy + 3
def f(x, y):
    return 4*x*y+3

# Ohjelman sisäiset parametrit
N = 1000
b = 20
s = 0.1

Ohjelma

Ohjelman varsinainen koodi ja logiikka on tässä osiossa.

In [23]:
# Tarvittavat kirjastot
import numpy as np
import matplotlib.pyplot as plt

# Arvotaan satunnaislukuja normaalijakaumasta
x = np.random.normal(mean_x, sd_x, N)
y = np.random.normal(mean_y, sd_y, N)

# Lasketaan funktion arvot N kertaa ja tallennetaan ne listaan
temp = []
for i in range(0, N - 1):
    temp.append(f(x[i], y[i]))

samples = np.array(temp)

# Lasketaan tulosjakauman tunnusluvut
mean = f(mean_x, mean_y) 
sd = 4 * np.sqrt(sd_x**2 * mean_y**2 + sd_y**2 * mean_x**2)

# Normaalijakauman tiheysfunktio
def pdf(x):
    return 1/(sd * np.sqrt(2 * np.pi)) * np.exp(- (x - mean)**2 / (2 * sd**2))

# Histogrammi
hist, bins = np.histogram(samples, bins=b) # muodostetaan lukumäärähistogrammi
width = bins[1] - bins[0]
area = np.sum(hist * width)
frequency = hist / area # lasketaan frekvenssit

# Lasketaan jakauman todelliset tunnusluvut
trueMean = np.mean(frequency) * area
trueSd = np.var(frequency) * area

deltaMean = round(abs(mean - trueMean) / mean, 2)
deltaSd = round(abs(sd - trueSd) / sd, 2)

Tulokset

Tulostetaan ohjelman laskemat tulokset. Histogrammin lisäksi ohjelma tulostaa muutamia arvottuja x-, y- ja laskettuja f(x,y)-arvoja. Lisäksi ohjelma tarkistaa, että histogrammi on normitettu oikein eli että pylväiden yhteispinta-ala on 1 (tämä ehto täytyy toteutua, jotta tiheysfunktion kuvaaja täsmää histogrammin kanssa). Lopuksi ohjelma vertaa lasketun ja todellisen jakauman tunnuslukuja.

In [24]:
# Yhteenveto

# Tulostetaan muutamia (x, y)-arvoja
print("Arvotaan", N, "normaalijakautunutta (x,y)-arvoa:")
print("x:     ", x[0:5], "...")
print("y:     ", y[0:5], "...")
print("f(x,y):", samples[0:5], "...")

# Piirretään histogrammi
center = (bins[:-1] + bins[1:]) / 2
plt.bar(center, frequency, width=width)
x = np.arange(bins[0], bins[-1], s)
plt.plot(x, pdf(x), color='r')
plt.title("Simuloinnin histogrammi")
plt.xlabel("f(x, y)")
plt.ylabel("Frekvenssi")
plt.show()

# Varmistetaan, että histogrammi on normitettu oikein 
# eli tarkistetaan, että pylväiden ala ≈ 1
sum = 0.0
for i in range(0, len(bins) - 1):
    sum += frequency[i] * width
sum = round(sum, 1)
print("Pylväiden pinta-ala:", sum, "OK" if sum - 1 < 0.1 else "- ei täsmää")

# Verrataan laskettua ja todellista jakaumaa
print("Laskettu keskikohta:", round(mean, 2), 
      ", todellinen keskikohta:", round(trueMean, 2), 
      ", ero:", deltaMean * 100, "%")
print("Laskettu leveys:", round(sd, 2), 
      ", todellinen leveys:", round(trueSd, 2), 
      ", ero:", deltaSd * 100, "%")
Arvotaan 1000 normaalijakautunutta (x,y)-arvoa:
x:      [1.57373097 1.42278662 1.55284284 1.95349008 1.7451489 ] ...
y:      [2.50299982 2.73210198 2.61514731 2.40050471 2.85259921] ...
f(x,y): [18.75619334 18.54879249 19.24365112 21.75744862 22.91284154] ...
Pylväiden pinta-ala: 1.0 OK
Laskettu keskikohta: 20.68 , todellinen keskikohta: 49.95 , ero: 142.0 %
Laskettu leveys: 2.32 , todellinen leveys: 2.59 , ero: 12.0 %
In [ ]: