1# -*- coding: utf-8 -*- 

2""" 

3@file 

4@brief Functions for :ref:`l-exemple_optim_alea`. 

5""" 

6import math 

7import random 

8 

9 

10def factorielle(x): 

11 """ 

12 Calcule :math:`x!` de façon récursive. 

13 """ 

14 if x == 0: 

15 return 1 

16 else: 

17 return x * factorielle(x - 1) 

18 

19 

20def profit(N, X, p, q, s): 

21 """ 

22 Calcule le profit. 

23 

24 @param N nombre de poulets vendus 

25 @param X nombre de poulets achetés 

26 @param p prix d'achat 

27 @param q prix de vente 

28 @param s prix soldé 

29 @return profit 

30 """ 

31 if X <= N: 

32 return X * (q - p) 

33 else: 

34 return X * (s - p) + N * (q - s) 

35 

36 

37def proba_poisson(lx, i): 

38 """ 

39 Calcule la probabilité :math:`\\pr{X=i}`` 

40 lorsque :math:`X` suit une loi de Poisson de paramètre 

41 :math:`\\lambda`. 

42 """ 

43 return math.exp(-lx) * (lx ** i) / factorielle(i) 

44 

45 

46def esperance(X, p, q, s, lx): 

47 """ 

48 Espérance du profit en faisant varier 

49 le nombre de poulet vendus. 

50 

51 @param X nombre de poulets achetés 

52 @param p prix d'achat 

53 @param q prix de vente 

54 @param s prix soldé 

55 @param lx paramètre :math:`\\lambda` 

56 @return espérance du profit 

57 """ 

58 res = 0.0 

59 for i in range(0, lx * 2): 

60 res += profit(float(i), X, p, q, s) * proba_poisson(lx, i) 

61 return res 

62 

63 

64def maximum(p, q, s, lx): 

65 """ 

66 Calcule les espérances de profit pour différents nombres 

67 de poulets achetés. 

68 

69 @param p prix d'achat 

70 @param q prix de vente 

71 @param s prix soldé 

72 @param lx paramètre :math:`\\lambda` 

73 @return liste ``(X, profit)`` 

74 """ 

75 res = [] 

76 for X in range(0, 2 * lx): 

77 r = esperance(X, p, q, s, lx) 

78 res.append((X, r)) 

79 return res 

80 

81 

82def find_maximum(res): 

83 """ 

84 Trouver le couple (nombre de poulets achetés, profit) 

85 lorsque le profit est maximum. 

86 

87 @param res résultat de la fonction :func:`maximum <mlstatpy.garden.poulet.maximum>` 

88 @return ``(X, profit)`` maximum 

89 """ 

90 m = (0, 0) 

91 for r in res: 

92 if r[1] > m[1]: 

93 m = r 

94 return m 

95 

96 

97def exponentielle(lx): 

98 """ 

99 Simule une loi exponentielle de paramètre :math:`\\lambda`. 

100 """ 

101 u = random.random() 

102 return - 1.0 / lx * math.log(1.0 - u) 

103 

104 

105def poisson(lx): 

106 """ 

107 Simule une loi de Poisson de paramètre :math:`\\lambda`. 

108 """ 

109 s = 0 

110 i = 0 

111 while s <= 1: 

112 s += exponentielle(lx) 

113 i += 1 

114 return i - 1 

115 

116 

117def poisson_melange(params, coef): 

118 """ 

119 Simule une variable selon un mélange de loi de Poisson. 

120 

121 @param params liste de paramètre :math:`\\lambda` 

122 @param coef ``coef[i]`` coefficient associé à la loi de paramètre ``params[i]`` 

123 @return valeur simulée 

124 """ 

125 s = 0 

126 for i, pa in enumerate(params): 

127 p = poisson(pa) 

128 s += p * coef[i] 

129 return s 

130 

131 

132def histogramme_poisson_melange(params, coef, n=100000): 

133 """ 

134 Calcule un histogramme d'un mélange de loi de Poisson. 

135 

136 @param params liste de paramètre :math:`\\lambda` 

137 @param coef ``coef[i]`` coefficient associé à la loi de paramètre ``params[i]`` 

138 @return histogramme 

139 """ 

140 h = [0.0 for i in range(0, 4 * max(params))] 

141 for i in range(0, n): 

142 x = poisson_melange(params, coef) 

143 if x < len(h): 

144 h[x] += 1 

145 s = sum(h) 

146 for i in range(0, len(h)): 

147 h[i] = float(h[i]) / s 

148 return h 

149 

150 

151def f_proba_poisson_melange(): 

152 """ 

153 Wraps function *proba_poisson_melange* to avoid 

154 global variable. 

155 """ 

156 

157 proba_poisson_melange_tableau = [] 

158 

159 def local_proba_poisson_melange(params, coef, i): 

160 """ 

161 Calcule la probabilité :math:`\\pr{X=i}`` 

162 lorsque :math:`X` suit un mélange de lois. 

163 

164 @param params liste de paramètre :math:`\\lambda` 

165 @param coef ``coef[i]`` coefficient associé à la loi de paramètre ``params[i]`` 

166 @return valeur 

167 """ 

168 if len(proba_poisson_melange_tableau) == 0: 

169 proba_poisson_melange_tableau.extend( 

170 histogramme_poisson_melange(params, coef)) 

171 if i >= len(proba_poisson_melange_tableau): 

172 return 0.0 # pragma: no cover 

173 return proba_poisson_melange_tableau[i] 

174 

175 return local_proba_poisson_melange 

176 

177 

178proba_poisson_melange = f_proba_poisson_melange()