← Back to branch summary

~ubuntu-branches/ubuntu/karmic/model-builder/karmic

« back to all changes in this revision

Viewing changes to model_builder/Bayes/lhs.py

  • Committer: Bazaar Package Importer
  • Author(s): Varun Hiremath
  • Date: 2007-04-10 17:05:04 UTC
  • Revision ID: james.westby@ubuntu.com-20070410170504-y884ntvt656218me
Tags: upstream-0.4.0
Import upstream version 0.4.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
model_builder/Bayes/lhs.py
 
1
#-----------------------------------------------------------------------------
 
2
# Name:        lhs.py
 
3
# Purpose:     Implements the Latin Hypercube Sampling technique as described by
 
4
#              Iman and Conover, 1982, including correlation control both for no 
 
5
#              correlation or for a specified correlation matrix for the sampled 
 
6
#              parameters.
 
7
#
 
8
# Author:      Antonio Guilherme Pacheco, Flavio Codeco Coelho
 
9
#
 
10
# Created:     2006/08/27
 
11
# RCS-ID:      $Id: lhs.py $
 
12
# Copyright:   (c) 2004
 
13
# Licence:     gpl
 
14
# New field:   Whatever
 
15
#-----------------------------------------------------------------------------
 
16
#!/usr/bin/python
 
17
"""
 
18
Implements the Latin Hypercube Sampling technique as described
 
19
by Iman and Conover, 1982, including correlation control
 
20
both for no correlation or for a specified correlation matrix 
 
21
for the sampled parameters
 
22
"""
 
23
from pylab import plot, figure,hist,show
 
24
#import random
 
25
import scipy.stats as stats
 
26
import numpy
 
27
from numpy.linalg import cholesky,inv
 
28
 
 
29
#import psyco
 
30
#psyco.full()
 
31
valid_dists = ['Normal','Triangular','Uniform']
 
32
def rank_restr(nparms=4, Iter=100, noCorrRestr=True, matCorr=None):
 
33
    """
 
34
    Does the correlation control.
 
35
    """
 
36
    x=[numpy.arange(1,Iter+1)for i in xrange(nparms)]
 
37
    if noCorrRestr:
 
38
        
 
39
        for i in xrange(nparms):
 
40
            numpy.random.shuffle(x[i])
 
41
            #x.append(numpy.array(random.sample(xrange(1,Iter+1), Iter)))
 
42
        return x
 
43
    else:
 
44
        if matCorr==None:
 
45
            C=numpy.core.numeric.identity(nparms)
 
46
        else:
 
47
            C=numpy.matrix(matCorr)
 
48
        s0=numpy.arange(1.,Iter+1)/(Iter+1.)
 
49
        s=stats.norm().ppf(s0)
 
50
        s1=[]
 
51
        for i in xrange(nparms):
 
52
            numpy.random.shuffle(s)
 
53
            s1.append(s.copy())
 
54
        S=numpy.matrix(s1)
 
55
        #print S,S.shape
 
56
        E=numpy.corrcoef(S)
 
57
        P=cholesky(C)
 
58
        Q=cholesky(E)
 
59
        Final=S.transpose()*inv(Q).transpose()*P.transpose()
 
60
        for i in xrange(nparms):
 
61
            x.append(stats.stats.rankdata(Final.transpose()[i,]))
 
62
        return x
 
63
 
 
64
def sample_cum(Iter=100):
 
65
    """
 
66
    Calculates samples on the quantile axis
 
67
    One sample per quantile
 
68
    """
 
69
    perc = numpy.arange(0,1,1./Iter)
 
70
    smp = [stats.uniform(i,1./Iter).rvs()[0] for i in perc]
 
71
    return smp
 
72
 
 
73
 
 
74
 
 
75
def sample_dist(cum, dist='Normal', parms=[0,1]):
 
76
    """
 
77
    Return a sample from from the given distribution 
 
78
    by inverting the cumulative prob. density function
 
79
    """
 
80
    if dist == 'Normal':
 
81
        if len(parms) == 2:
 
82
            n = stats.norm(parms[0],1./parms[1])
 
83
            d = n.ppf(cum)
 
84
    elif dist=='Triangular':
 
85
        if len(parms) ==3 and parms[0]<=parms[1]<=parms[2]:
 
86
            loc=parms[0]
 
87
            scale=parms[2]-parms[0]
 
88
            t = stats.triang((float(parms[1])-loc)/scale,loc=loc,scale=scale)
 
89
            d = t.ppf(cum)
 
90
    elif dist == 'Uniform':
 
91
        loc = parms[0]
 
92
        scale = parms[1]-parms[0]
 
93
        d = stats.uniform(loc,scale).ppf(cum)
 
94
    else:
 
95
        print '%s is an unsupported distribution!'%dist
 
96
        
 
97
    return d
 
98
 
 
99
 
 
100
 
 
101
def lhs(Pars, dists, parms, Iter=100, noCorrRestr=True, matCorr=None):
 
102
    """
 
103
    Returns tuple of vectors
 
104
    Pars: list of strings with parameter names
 
105
    dists: List of strings with distribution names
 
106
    Iter: sample size
 
107
    defaults to no correlation between sampled parameters
 
108
    """
 
109
    ParsList=[]
 
110
    if len(Pars)==len(dists):
 
111
        indexes=rank_restr(nparms=len(dists), Iter=Iter, noCorrRestr=noCorrRestr, matCorr=matCorr)
 
112
        for i in xrange(len(dists)):
 
113
            v=sample_dist(sample_cum(Iter), dist=dists[i], parms=parms[i])
 
114
            index=map(int,indexes[i]-1)
 
115
            ParsList.append(v[index])
 
116
    return ParsList
 
117
            
 
118
    
 
119
 
 
120
if __name__=='__main__':
 
121
    #c=lhs(['Par1', 'Par2', 'Par3'],['Normal','Triangular','Uniform'], [[0,1], [1,5,8], [1,2]],100000)
 
122
    c=lhs(['Par1', 'Par2', 'Par3'],['Normal','Triangular','Uniform'], [[0,1], [1,5,8], [1,2]],10000, noCorrRestr=True)
 
123
  #  m=[[1,.3595,-.2822],[.3595,1,-.1024],[-.2822,-.1024,1]]
 
124
  #  c=lhs(['Par1', 'Par2', 'Par3'],['Normal','Triangular','Uniform'], [[0,1], [1,5,8], [1,2]],100000, noCorrRestr=False, matCorr=m)
 
125
    print stats.spearmanr(c[0],c[1])
 
126
    print stats.spearmanr(c[0],c[2])
 
127
    print stats.spearmanr(c[1],c[2])
 
128
    
 
129
##    print c
 
130
##    print 'done!'
 
131
    #a=sample_cum(10000)
 
132
    #b=sample_dist(a, dist='Triangular', parms=[1,5,8])
 
133
    
 
134
    #plot(b,a, 'bo')
 
135
##    print c[0].shape,c[1].shape,c[2].shape
 
136
##    print c[0][0], type(c[0][0])
 
137
    hist(c[0], bins=30)
 
138
    figure()
 
139
    hist(c[1], bins=30)
 
140
    figure()
 
141
    hist(c[2], bins=30)
 
142
    show()
 
143