truncnorm() in python

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
asahdkhaled
User
Beiträge: 29
Registriert: Samstag 28. Oktober 2017, 22:07

Hallo,
ich habe folgendes Problem:
Ich möchte für ein Intervall, mit der unteren Grenze a und der oberen Grenze b, µ (in der Mitte des Intervalls) und der standardabweichung Daten generieren, und zwar per "truncated Gaussian distribution".

Das was ich bisher gefunden habe ist:

Code: Alles auswählen

myclip_a = 0
myclip_b = 1
my_mean = 0.5
my_std = 0.3

a, b = (myclip_a - my_mean) / my_std, (myclip_b - my_mean) / my_std
x_range = np.linspace(-1,2,1000)
plt.plot(x_range, truncnorm.pdf(x_range, a, b, loc = my_mean, scale = my_std))
Das beantwortet aber nur beschränkt mein Problem.
Wie man sieht, ist es eher eine normal Verteilung die an beiden Enden des Intervalls abgeschnitten ist.
Ich möchte aber innerhalt dieses kleinen Intervalls, eine eigene Normalverteilung Programmieren.
Wie geht das?

truncnorm() ist wohl meine Wahl#1, allerdings blicke ich noch nicht richtig durch die Parameter durch.
truncnorm(a,b,µ,sd) wäre für mich sinnvoller, kann mir da jemand auf die Sprünge helfen?
narpfel
User
Beiträge: 643
Registriert: Freitag 20. Oktober 2017, 16:10

Moin,

was genau meinst du mit „innerhalt dieses kleinen Intervalls, eine eigene Normalverteilung“? `truncnorm.pdf` ist eine Normalverteilung in dem angegebenen Intervall, normiert dass ihr Integral 1 ist.

Falls dich die Umrechnung
asahdkhaled hat geschrieben:

Code: Alles auswählen

a, b = (myclip_a - my_mean) / my_std, (myclip_b - my_mean) / my_std
stört, kannst du ja einfach eine Funktion benutzen, die das für dich erledigt:

Code: Alles auswählen

def truncated_gaussian_pdf(x, lower, upper, mu, sigma):
    a = (lower - mu) / sigma
    b = (upper - mu) / sigma
    return truncnorm.pdf(x, a, b, loc=mu, scale=sigma)

# Test
np.allclose(
    truncated_gaussian_pdf(x, lower, upper, mu, sigma),
    truncnorm.pdf(x, a, b, mu, sigma)
)
`truncnorm.pdf` ist ungefähr so implementiert:

Code: Alles auswählen

def truncated_gaussian_pdf2(x, lower, upper, mu, sigma):
    y = norm.pdf(x, mu, sigma)
    y[x < lower] = 0
    y[x > upper] = 0
    y /= np.diff(norm.cdf([lower, upper], mu, sigma))
    return y

# Test
np.allclose(
    truncnorm.pdf(x, a, b, mu, sigma),
    truncated_gaussian_pdf2(x, lower, upper, mu, sigma)
)
Hilft das bei der Beantwortung deiner Frage?
asahdkhaled
User
Beiträge: 29
Registriert: Samstag 28. Oktober 2017, 22:07

Ich meine, dass ich sagen wir mal 3 Intervalle habe:
1. [0,0.4]
2.[0.4,0,8]
3.[0,8,1]

Innerhalb dieser 3 Intervalle möchte ich Datensätze samplen, die truncated gaussian distribution zur Grundlage haben.

Bsp für Intervall 1:

truncnorm(0,0.4,µ(von 1), sd (von1).

Wo ist der unterschied zwischen truncnorm() und truncnorm().pdf??
Außerdem verstehe ich nicht, was x_range bringen soll? Und was ist der Unterschied zwischen a, b und den lower und upper Grenzen?
narpfel
User
Beiträge: 643
Registriert: Freitag 20. Oktober 2017, 16:10

Ok, du möchtest also gar nicht mit der PDF arbeiten, sondern Zufallszahlen erzeugen, die truncnorm-verteilt sind? Dafür gibt es `truncnorm.rvs`.

Zum Unterschied von `truncnorm` `truncnorm.pdf`: Das ist in der Dokumentation zu `scipy.stats.rv_continuous` beschrieben (Punkt „Frozen Distributions“).

`x_range` sind die Werte, an denen du deine PDF auswertest.

Der Unterschied zwischen `a`, `b` und `lower` bzw. `upper` (bei dir im Code `myclip_{a, b}`) wird deutlich, wenn du dir klar machst, wie die truncated PDF berechnet wird: `truncnorm` verwendet intern die PDF der Normalverteilung. Für (μ, σ) ≠ (0, 1) wird einfach eine Streckung und Verschiebung der x-Achse durchgeführt: x → (x - μ) / σ. Dabei werden die Intervallgrenzen vor der Koordinatentransformation an die PDF angelegt und werden somit mittransformiert. Man muss also die gewünschten Grenzen so angeben, dass sie nach der Transformation stimmen. Genau das macht die Umrechnung. Deutlich wird das am ehesten, wenn man sich mal ein paar verschiedene PDFs plottet und sich die Auswirkungen der einzelnen Parameter ansieht.
asahdkhaled
User
Beiträge: 29
Registriert: Samstag 28. Oktober 2017, 22:07

Ich habe nochmal dazu eine Frage:

Die truncnorm() Function erwartet ja den Mittelwert und die Standardabweichung als Parameter für das Intervall. Wie kann ich das ausrechnen für dasspezielle Intervall in Python?Die Funktion numpy.mean() scheint nicht zu klappen, ich bekomme komische Ergebnisse, die Ergebnisse die truncnorm() zurück gibt, sind zu großen Teilen außerhalb des Intervalls. Das ist besonders bei Intervallen die sehr klein sind, wo vll die richtige Berechnung des Mittelwerts und Standardabweichung sehr wichtig ist.
Für breitere Intervalle klappt es eg ganz gut. Gibt es ein Limit wie klein ein Intervall sein darf?

Beispiel:

[0,12;0,17]--> Value 0,0937818650369 (out of range)*

Ich copy/paste auch mal meinen Forenbeitrag aus Stackoverflow, ich hoffe es ist ok, dass der Beitrag in english ist;


What I want to do is: I have an Intervall, sample me one Value, which is between the boarder of that intervall and simple it in a way of truncated normal distribution. I have an extra column and it should write the value I gain by sampling in another column. For example: Intervall [0.2;0.6] --> sample value 0.343433 I think I found a solution:

truncnorm().stats()

But I don't know why, but for the parameters I give the

truncnorm()

function, almost 50% of the values I gain are outside the boarders. What am I doing wrong?

Here is the code (a short part of the code)

Code: Alles auswählen

      convert_cat=(name_convert_column,name_convert_column,_tabelle,name_convert_column,_tabelle,_tabelle,name_convert_column)
    drop_view=(name_convert_column)
    calculate=(name_convert_column,name_convert_column,name_convert_column,name_convert_column,name_convert_column,_tabelle,name_convert_column,name_convert_column)
    cur.execute("CREATE VIEW convert_cat_%s (quotient, %s, rnum) AS SELECT (COUNT(*)/(SELECT COUNT(*) FROM %s ) ) as quotient, %s, row_number() over ( order by (COUNT(*)/(SELECT COUNT(*) FROM %s ) ) desc ) as rnum FROM     %s  GROUP BY %s ORDER BY quotient desc" %convert_cat)
    cur.execute("Select b.ID,a.unten,a.oben, a.mean, a.sd FROM( SELECT t3.RNUM, t3.%s, lag(t3.com_Pr,1,0) OVER (order by rnum asc) as unten , t3.com_PR as oben, ((t3.com_PR +(lag(t3.com_Pr,1,0) OVER (order by rnum asc)))/2) as MEAN, ((t3.com_PR-(lag(t3.com_Pr,1,0) OVER (order by rnum asc)))/6) AS SD FROM( SELECT t1.rnum, t1.%s , SUM(t2.quotient) as com_Pr FROM CONVERT_CAT_%s t1 INNER JOIN CONVERT_CAT_%s t2 ON t1.rnum >= t2.rnum group by t1.rnum, t1.%s, t1.quotient ORDER BY RNUM asc ) t3) a INNER JOIN %s b ON b.%s = a.%s order by ID asc" %calculate)
    _content_category = cur.fetchall()
    add_category_number_column = (_tabelle, name_convert_column)
    cur.execute("ALTER TABLE %s ADD %s_category NUMBER(15,14)" % add_category_number_column)
    x=0
    for ID in _content_category:
        id = _content_category[0]
        id_category = [j[0] for j in _content_category]
        unten_category = [j[1] for j in _content_category]
        oben_category = [j[2] for j in _content_category]
        #mean_category = [j[3] for j in _content_category]
        sd_category = [j[4] for j in _content_category]
        mean, var = truncnorm.stats(unten_category[x], oben_category[x], moments='mv')
       # sd = np.sqrt(var)
        X = get_truncated_normal(mean= mean, sd=sd_category[x], low=unten_category[x], upp=oben_category[x])
        update_cells_value = float(X.rvs(1))
        category = (_tabelle, name_convert_column,update_cells_value,id_category[x])
     cur.execute("UPDATE %s SET %s_category = %s WHERE ID=%s" % category)

        x += 1
I tried to calculate mean and sd in the sql query with

1) ((t3.com_PR +(lag(t3.com_Pr,1,0) OVER (order by rnum asc)))/2) as MEAN
2) ((t3.com_PR-(lag(t3.com_Pr,1,0) OVER (order by rnum asc)))/6) AS SD

and with truncnorm().stats() function. Seems that with the stats function, the result are getting even worse and the values are even more out of range then before...
narpfel
User
Beiträge: 643
Registriert: Freitag 20. Oktober 2017, 16:10

Zuerst mal hast du ein ziemlich großes Problem mit deinem SQL-Code: Der ist anfällig für SQL-Injections.

Ich bin mir nicht sicher, dass ich verstehe, was du machen willst. Den Mittelwert und die Standardabweichung einer Verteilung muss man kennen, bevor man da irgendwas damit machen kann. Ohne diese beiden Werte haben die Zahlen, die am Ende herauskommen, überhaupt keine Bedeutung...

Was versuchst du denn eigentlich zu tun?

Ansonsten gilt das, was ich schon geschrieben habe: Wenn (μ, σ) ≠ (0, 1) ist, dann musst du deine Intervallgrenzen entsprechend umrechnen. Dann kommen da auch keine Werte außerhalb des Intervalls raus.


Link zum Crosspost.
asahdkhaled
User
Beiträge: 29
Registriert: Samstag 28. Oktober 2017, 22:07

Danke dir erstmal für deine fleißige Hilfe!
Das mit den Grenzen umrechnen habe ich wie von dir vorgeschlagen mit

Code: Alles auswählen

def get_truncated_normal(mean, sd, low, upp):
    return truncnorm((low - mean) / sd, (upp - mean) / sd, loc=mean, scale=sd)
  

gelöst. Ist da viellicht noch irgendetwas zu beachten?
Das Ziel soll es sein, dass ich eine x-Achse habe, in dieser x-Achse gibt es verschiedene Intervalle.
Inerhalb dieser Intervalle will ich zufällig eine Zahl generieren (Bspw. Intervall [1;3] und dann soll er bspw. 2 auswerfen.
Bei einem weiteren Versuch dann 1.2 oder 2.9 usw usw.

Das Problem ist dadurch, dass die Intervalle ja jeweils eine eigene Standardabweichung/Mittelwert haben und ich den ja anscheinend erstmal bestimmten muss, bevor ich die truncnorm() Funktion anwende. Dazu muss ich ja dann zunächst erstmal mit den Grenzen und dem wissen das es sich um eine truncated normal distribution handelt, Mittelwert und SD ausrechnen können. Kann mir dabei jemand helfen? Die bisherigen Vorschläge gehen ja eher davon aus, dass mean und sd schon bekannt ist...
narpfel
User
Beiträge: 643
Registriert: Freitag 20. Oktober 2017, 16:10

Ich wiederhole meine Frage: Was willst du eigentlich tun? Was bedeutet die Zahl, die du generieren willst? Was ist das für ein Intervall? Willst du vielleicht eine gleichverteilte Zahl (im Gegensatz zu einer (trunc-) normalverteilten)? Ist das hier eine Instanz eines XY-Problems?

Wenn du einfach mal ein paar PDFs von Truncnormal-Verteilungen mit verschiedenen (beliebigen) μs und σs plottest, hilft dir das bei deinem Problem?
asahdkhaled
User
Beiträge: 29
Registriert: Samstag 28. Oktober 2017, 22:07

Also was ich vorhabe:

1. Ich habe eine Spalte in meiner DB die categorical verteilt ist (Beispiel Geschlecht: Männlich, Weiblich)
2. Dann berechne ich die kummulierte Wahrscheinlichkeitsverteilung und Teile dann eine Kategorie einem bestimmten Intervall zu.
Bspw: von 100 Datensätzen 70 weiblich, 30 männlich
p(Männlich)=0,3
p(weiblichnich)=0,7
[0;0,7] Ist das Intervall für weiblich
und [0,7;1] für männlich.

Jetzt möchte ich mir für die jeweiligen Intervalle eine Zufallszahl die normalverteilt ist generieren lassen, die bspw. für weiblich einen Wert zwischen 0 und 0,7 zurück gibt. Dabei komplett random. Da ich mehrere Spalten habe mit denne ich das machen will und die teilweise um die 20 Kategorien haben, habe ich zuvor also noch nicht die SD und MEAN.

Wenn die Intervalle kleiner werden, sprich mehr Kategorien exisitieren, in die dann also das Intervall [0;19 aufgteilt wird, kommt es zu Problem und es gibt Datensätze die außerhalb des Intervalls liegen.

ich will also nicht eine Normalverteilung einfach nur abschneiden und damit eine an den Seiten abgeschnittene Normalverteilung haben, sondern tatsächlich für jedes einzelne Intervall (welches eine Kategorie wiederspiegelt) eine eigene Normalvertielungsfunktion.

Konttst du etwas damit anfangen? Da nach wie vor Fehler auftauchen, auch wenn ich die Werte umrechne, gehe ich momentan also davon aus, dass die Mittelwert bzw. SD Berechnung für das jeweilige Intervall nicht hinhaut
__deets__
User
Beiträge: 14480
Registriert: Mittwoch 14. Oktober 2015, 14:29

Ich bin nur so milde bewandert in Statistik, aber ist deine Modellierung nicht Unfug? Eine Normalverteilung impliziert doch, das du immer einen Wert bekommen kannst, der mit einer gewissen Wahrscheinlichkeit auch irgendwo weit außerhalb liegt.

Machen wir ein Gedanken Experiment: du setzt u auf den Mittelwert deines Intervalls, und s auf die breite/2. Dann liegen ~60% deiner Werte im Intervall. Die restlichen 40 % musst du clampen. Jetzt kannst du s natürlich beliebig verkleinern, und irgendwann hast du nur eine geringe Chance, das der Wert außerhalb des Intervalls liegt. Passieren kann es aber immer noch, du musst also clampen (bzw truncaten).

Das mean für das Mittel des Intervalls ist ja nun trivial auszurechnen. Nur bei der std wirst du eine Vorgabe machen müssen, relativ zur breite. Das ist aber deine eigene Entscheidung, da gibt es nichts richtiges oder falsches.
asahdkhaled
User
Beiträge: 29
Registriert: Samstag 28. Oktober 2017, 22:07

Naja, ich bin auch kein Experte, aber wenn ich doch sage ich mal die Grenzen setze und innerhalb dieser Grenze eine Normalverteilung haben möchte, wo der minimale Wert den die truncnorm() zurück gibt, auf der Intervall Grenze ist, muss das ja irgendwie möglich sein?!
Also der mean ist klar, für die Standardabweichung habe ich auch eine Vorgabe, (obere_grenze-untere_Grenze)/6.

Aber ich finde es halt komisch, dass man, obowhl man Grenzwerte angibt, trotzdem Daten außerhalb dieser Grenze ausgespuckt werden :K
narpfel
User
Beiträge: 643
Registriert: Freitag 20. Oktober 2017, 16:10

Okay, das ist jetzt eine neue Information, dass du doch einen Wert für σ gegeben hast. Ich kann allerdings immer noch nicht nachvollziehen, ob das ganze so Sinn macht...

Bei mir kommen keine Werte außerhalb des Intervall vor:

Code: Alles auswählen

lower = 0.12
upper = 0.13
μ = (lower + upper) / 2
σ = (upper - lower) / 6
a = (lower - μ) / σ
b = (upper - μ) / σ

samples = truncnorm.rvs(a, b, μ, σ, 1_000_000)
np.any((samples < lower) | (samples > upper))
# → False
asahdkhaled
User
Beiträge: 29
Registriert: Samstag 28. Oktober 2017, 22:07

Ich kann es mir zwar gerade auch nicht erklären, aber bei mir funktioniert es jetzt auch :) (Glaube ich zumindest :D )
__deets__
User
Beiträge: 14480
Registriert: Mittwoch 14. Oktober 2015, 14:29

Ich vermute mal du hast die Umrechnung des a/b Intervalls vergessen. Das was narpfel in Zeile 5 und 6 tut.
narpfel
User
Beiträge: 643
Registriert: Freitag 20. Oktober 2017, 16:10

asahdkhaled hat geschrieben:

Code: Alles auswählen

def get_truncated_normal(mean, sd, low, upp):
    return truncnorm((low - mean) / sd, (upp - mean) / sd, loc=mean, scale=sd)
Naja, hier ist die Umrechnung noch da. Parameter vertauscht wäre auch eine Möglichkeit. :wink:
__deets__
User
Beiträge: 14480
Registriert: Mittwoch 14. Oktober 2015, 14:29

Aber er füttert das Ding doch danach. Wenn er das nicht mit ebenfalls adaptierten Intervallgrenzen tut, kommt da auch Unsinn raus.
asahdkhaled
User
Beiträge: 29
Registriert: Samstag 28. Oktober 2017, 22:07

Kommando zurück, komischerweise kommen wieder falsche Werte raus :D Bin also genauso schlau wie vorher
__deets__
User
Beiträge: 14480
Registriert: Mittwoch 14. Oktober 2015, 14:29

Wie fütterst du das denn? Zeig mal ein bisschen mehr statt nur die Umrechnungsfunktion.
narpfel
User
Beiträge: 643
Registriert: Freitag 20. Oktober 2017, 16:10

@__deets__: Das von `truncnorm` zurückgegebene Objekt hat alle Parameter der Verteilung fest eingebaut. Wenn man darauf später z. B. `rvs` aufruft, dann muss man die Parameter nicht nochmal angeben.

@asahdkhaled: Mehr Code wäre tatsächlich hilfreich.
asahdkhaled
User
Beiträge: 29
Registriert: Samstag 28. Oktober 2017, 22:07

Ich blicke nicht mehr durch. Eben waren noch alle Werte korrekt. Jetzt nochmal durchlaufen lassen, jetzt sind wieder ca. 50% falsch....
Der Code ist fast der gleiche wie vorher, alles was zuvor passiert hat damit nichts zu tun...)

Code: Alles auswählen

for datatypes in datentypen_columns:
    name_convert_column = results_as_list_new[i]
    if datatypes =='VARCHAR2':
        convert_cat=(name_convert_column,name_convert_column,_tabelle,name_convert_column,_tabelle,_tabelle,name_convert_column)
        drop_view=(name_convert_column)
        calculate=(name_convert_column,name_convert_column,name_convert_column,name_convert_column,name_convert_column,_tabelle,name_convert_column,name_convert_column)
        cur.execute("CREATE VIEW convert_cat_%s (quotient, %s, rnum) AS SELECT (COUNT(*)/(SELECT COUNT(*) FROM %s ) ) as quotient, %s, row_number() over ( order by (COUNT(*)/(SELECT COUNT(*) FROM %s ) ) desc ) as rnum FROM     %s  GROUP BY %s ORDER BY quotient desc" %convert_cat)
        cur.execute("Select b.ID,a.unten,a.oben FROM( SELECT t3.RNUM, t3.%s, lag(t3.com_Pr,1,0) OVER (order by rnum asc) as unten , t3.com_PR as oben FROM( SELECT t1.rnum, t1.%s , SUM(t2.quotient) as com_Pr FROM CONVERT_CAT_%s t1 INNER JOIN CONVERT_CAT_%s t2 ON t1.rnum >= t2.rnum group by t1.rnum, t1.%s, t1.quotient ORDER BY RNUM asc ) t3) a INNER JOIN %s b ON b.%s = a.%s order by ID asc" %calculate)
        _content_category = cur.fetchall()
        add_category_number_column = (_tabelle, name_convert_column)
        cur.execute("ALTER TABLE %s ADD %s_category NUMBER(15,14)" % add_category_number_column)
        x=0
        for ID in _content_category:
            id = _content_category[0]
            id_category = [j[0] for j in _content_category]
            lower = [j[1] for j in _content_category]
            upper = [j[2] for j in _content_category]
            mean_category = ((lower[x] + upper[x]) / 2)
            sd_category = (upper[x] - lower[x]) / 6
            a = (lower[x] - mean_category) / sd_category
            b = (upper[x] - mean_category) / sd_category
            samples = float(truncnorm.rvs(a, b, mean_category, sd_category, 1))
            #X = get_truncated_normal(mean= mean_category[x], sd=sd_category[x], low=unten_category[x], upp=oben_category[x])
            #update_cells_value = float(X.rvs(1))
            category = (_tabelle, name_convert_column,samples,id_category[x])
            cur.execute("UPDATE %s SET %s_category = %s WHERE ID=%s" % category)
            x += 1
        drop_view = (name_convert_column)
        #cur.execute("DROP VIEW convert_cat_%s" % drop_view)

Antworten