---

Hej!

Jeg har r1,r2,...,rN samplede værdier som definere en eller
anden kontinuert fordeling P(r). Hmm. giver dette udsagn
nogen som helst mening?

Jeg vil nu gerne udregne dP/dr /P = dlog(P(r))/dr og jeg vil gerne
udregne det uden at bestemme P, men udelukkende baseret på mine
sampled værdier (pga. æstetisk tilfredsstillelse).

Motivationen er at hvis p(r) er en kendt sandsynelighedsfordeling
for afstanden mellem to partikler, så er A=-TS
= -T*kb log(P(r))+konst den fri energi svarende til en Boltzmann
vægt P(r)=Nexp(-A(r)/kbT) og F=-dA/dr er gennemsnitskraften
mellem partiklerne i en statistisk fysisk forstand (A er
potential of mean force).

Jeg er interesseret i at gøre det samme for en polymer, og
der finde gennemsnitskraften mellem monomere i en bestemt
afstand langs polymeren, fordi dette kan bruges til at udregne
elastiske egenskaber. Til laveste orden er fordelingen gaussisk,
og gennemsnitskraften er derfor proportional med afstanden, dvs.
en hooksk fjeder med en fjederkonstant, der afhænger af hvor
langt de to monomere er fra hinanden langs polymerkæden.
Det interessante er netop afvigelserne fra denne opførsel,
som afstanden mellem monomere bliver kort, fordi det giver
de elastiske egenskaber på mikroskopiske skala.

Anyways det er motivationen, problemet er matematisk nogle
ideer hvordan det kan gøres på en "let måde"?

Min idee er ganske forfærdelig besværlig, nemmeligt at udregne
n momenter af fordelingen. Det giver <r^k> k=1,..,n

Når jeg har momenterne så kan en maximum entropi ansatz
bruges til at definere sandsynelighedsfordelingen, der
har formen: P(r)= exp(l[0] + l[1]r + l[2]r² + ..)
hvor l[0],..,l[n] Lagrange multipliers, der bestemmes ved
maksimering af entropien S[P] = integral -kP(r)lnP(r) +
de n constraints, der defineres af de momenter jeg får fra
de samplede værdier.

Det er så let at finde
dlog(P(r))/dr = l1 + 2*l2r +.. + n*l[n]*r^(n-1)

Hvor alle lagrange multipliers er implicitte funktioner
af de samplede momenter. Problemet er at det ikke er let
at finde lagrangemultipliers med denne fremgangsmetode.

Fordelingen P(r) er formodeligt tæt på en gaussisk, så
jeg behøver en 5 momenter for at få afvigelsen med, og
så bliver det en "interessant" integration for at finde
moments. Alternativet er numerisk minimation med
constraints hvilket også lyder som en "interessant" opgave.

Et alternativ kunne være en ekspansion med hermit
polynomier?

Nogle der har en mere praktisk idee?

--
Mvh. Carsten Svaneborg
http://www.softwarepatenter.dk where you do not
want to go in the future!

---

Carsten Svaneborg <zqex@nowhere.on.the.net> writes:
> Min idee er ganske forfærdelig besværlig, nemmeligt at udregne
> n momenter af fordelingen. Det giver <r^k> k=1,..,n

Er det virkeligt så langsomt at udregne momenterne? Det kræver n
gangeoperationer per datapunkt r^3 = r * r^2 og r^4 = r^2 * r^2
etc.. Det virker ikke så slemt, men jeg har ikke prøvet det?

I princippet kan du sikkert spare lidt tid ved at lave en tabel over
funktionsværedierne af nogle ortogonale polynomier (eller andre
funktioner). Hvis du bruger blitz++ (til c++) kan du måske finde en
metode til at lave tabellen i compile time, men så vidt jeg har
forstået ledte du efter en smuk løsning, og rekursive templates er
ikke særligt smukke :)

Egentlig tror jeg ikke at det giver store fejl ved at inddele r-aksen
i små intervaller, og derefter bruge intervalhyppighederne til at
regne momenterne ud. Så vidt jeg har forstået er denne metode standard
hvis man vil regne Van Hove funktioner ud, og det en det hjørne af
MD-simulering som jeg lige har kigget på. Der må står vist noget om
det i Allen og Tildesley, men jeg kan ikke lige finde den. Om ikke
andet er det let at undersøge hvor store fejl det giver ved at regne
momenterne ud på to måder (med og uden intervaller).

Skriv hvis du finder på en smuk løsning :)

Niels

PS: Er det egentligt hurtigt at caste til en variabeltype med færre
cifre? (long double til en float el lign)?

--
Niels L Ellegaard http://dirac.ruc.dk/~gnalle/

---

Niels Langager Ellegaard wrote:
>> Min idee er ganske forfærdelig besværlig, nemmeligt at udregne
>> n momenter af fordelingen. Det giver <r^k> k=1,..,n
> Er det virkeligt så langsomt at udregne momenterne?

Nej, det er ikke noget problem. For at få en idee om errorbars kører
jeg allerede analysen som 10 separate jobs, og hver tager mindre end
5 timer, tilgengæld kræver de 600-700Mb hver.

Problemet er at for et givet sæt momenter at udregne de tilsvarende
lagrange multipliers, fordi det er et maksimeringsproblem med N
constraints. Og det er værdierne af lagrange multipliers jeg skal
bruges til at udregne det effektive force field.

> Skriv hvis du finder på en smuk løsning :)

Ikke særligt æstetisk, men til at begynde med vil jeg blot lave
et histogram, og så fitte et polynomium til log(P), det er så
let at udregne kraftfeltet blot ved at differentiere det fittede
polynomium.

--
Mvh. Carsten Svaneborg
http://www.softwarepatenter.dk where you do not
want to go in the future!

---

Carsten Svaneborg <zqex@nowhere.on.the.net> writes:
> Niels Langager Ellegaard wrote:
> Nej, det er ikke noget problem. For at få en idee om errorbars kører
> jeg allerede analysen som 10 separate jobs, og hver tager mindre end
> 5 timer, tilgengæld kræver de 600-700Mb hver.
>
> Problemet er at for et givet sæt momenter at udregne de tilsvarende
> lagrange multipliers, fordi det er et maksimeringsproblem med N
> constraints. Og det er værdierne af lagrange multipliers jeg skal
> bruges til at udregne det effektive force field.

Hmm, det havde jeg vist misforstået. Her er en ny ide, der sikkert
ikke virker, men den var skæg at finde på, og den er måske hurtig at
teste. Definer den karakteristiske funktion phi(omega) ved

phi(omega) = \sum_n <(r - <r>)^n> / n! * (i omega)^n

Så gælder det så vidt jeg husker at phi(omega) er den komplekst
Fouriertransformerede af sandsynlighedstætheden.

phi(omega) = int_{-\infty}^\infty exp(i omega r) p(r) dr

Så måske kan du finde p(r) og p'(r) ved hjælp af ifft. Du kan jo bare
nøjes med dine 5-10 momenter og 1000 værdier af s.

Skriv hvis det virker. Jeg er selv i tvivl, men jeg har ikke testet
metoden selv (for det er jo usenet :) )

Niels

Her står lidt om den kakateristiske funktion:

http://mathworld.wolfram.com/CharacteristicFunctionProbability.html

--
Niels L Ellegaard http://dirac.ruc.dk/~gnalle/

---

Niels L. Ellegaard wrote:
> Hmm, det havde jeg vist misforstået. Her er en ny ide, der sikkert
> ikke virker, men den var skæg at finde på, og den er måske hurtig at
> teste. Definer den karakteristiske funktion phi(omega) ved
> phi(omega) = \sum_n <(r - <r>)^n> / n! * (i omega)^n

Dit moment er lidt mystisk, ifht. ligningen jeg husker men jeg
forstår hvad du mener.

Faktisk er det endnu bedrer at bruge kumulanter, fordi de giver
udviklingen af ln(phi(omega)) istededet for phi(omega), og den
afledede kan så udregnes analytisk i termer af den trunkerede
serie af kumulanter og kraf vs. afstand kan findes direkte.
Desuden er fordelingen gaussisk for n stor nok, hvorfor
højreordens kumulanterne forsvinder.

> Skriv hvis det virker. Jeg er selv i tvivl, men jeg har ikke testet
> metoden selv (for det er jo usenet :) )
Nu skal jeg lige på ferie først ... ;*)

--
Mvh. Carsten Svaneborg
http://www.softwarepatenter.dk where you do not
want to go in the future!

---

Niels Langager Ellegaard wrote:

> PS: Er det egentligt hurtigt at caste til en variabeltype
> med færre cifre? (long double til en float el lign)?

Det er vist en af de ting, der afhænger voldsomt af ens
processor og oversætter. Jeg ville simpelthen prøve det af.

--
Jens Axel Søgaard