Simulations with Bernoulli models

5 Simulations with Bernoulli models

In this section, we perform simulations of sequences build a Bernoulli models, and compute the log-likelihoods on them.

import partition import random import sequence import lcompte import matrice import parti_simp import compte import lexique

The function that generates a Bernoulli model:

def gen_bern(lalpha): lpr=lcompte.Lproportion() for i in range(len(lalpha)): alpha=lalpha[i] p = compte.Proportion() s="|A %f\n|B %f\n"%(alpha,1-alpha) p.read_str(s) lpr[i]=p return lpr

We define some parameters

lalpha=[0.3,0.7] # parameters of the models lseq=10000 #length of the sequences nsegprob=100 # length of the computation nseg=15 # nb of simulated segments nbsampl=100 # nb of samples

and generate the Lexique.

ldesc=lpr.num() lldesc=len(ldesc) if lldesc<=1: print "Trop few models" exit lx=lexique.Lexique(Lprop=lpr)

The output file:

s=reduce(lambda x,y: x+"_"+y, map(str, lalpha)) fprob=open("bern_%s_n%d.llh"%(s,nseg),"w") for i in range(nsegprob): fprob.write("N%d\t"%(i+1)) fprob.write("\n")

And the simulation loop:

for i in range(nbsampl): p=partition.Partition() #generation of the partition p.build_random(lseq,nseg,ec=50) ns=random.choice(ldesc) for sg in p: sg.g_num([ns]) ns2=ns while ns2==ns: ns=random.choice(ldesc) s=sequence.Sequence() s.generate(lseq) for sg in p: s.read_prop(lpr[sg.num()[0]],deb=sg.deb(),fin=sg.fin()) ############### # Log-likelihoods lp=lx.log_likelihood(s,nsegprob) ############### # output for x in lp: fprob.write("%.12f\t"%(x)) fprob.write("\n") fprob.flush() fprob.close()