1pub type CDHResult<T> = Result<T,String>;
16
17
18#[derive(Debug)]
20pub struct Probability;
21
22#[derive(Debug,Clone)]
28pub struct RandomVariable{
29pub variable:Vec<f64>,
30pub probability: Vec<f64>
31}
32
33#[derive(Debug)]
35pub struct RandomVector<T>{
36pub variables: Vec<RandomVariable>,
37pub probability_matrix:JointProbability<T>}
38
39#[derive(Debug,PartialEq)]
41pub struct JointProbability<T>{
42pub matrix:Vec<Vec<T>>
43}
44
45
46pub fn check_probability()
51->CDHResult<()>{
52println!("+++++ probability.rs results +++++");
53let x = (0..10).map(|i| i as f64).collect::<Vec<f64>>();
54let px = x.iter().map(|_| 1f64/10f64)
55.collect::<Vec<f64>>();
56let rv = RandomVariable::new(x,px);
57println!("Random Variable Mean === {:?}"
58,rv.mean()?);
59println!("Random Variable Variance === {:?}"
60,rv.variance()?);
61println!("Random Variable Standard Deviation === {:?}"
62,rv.std_deviation()?);
63println!("Random Variable Moment Generating func(5) === {:?}"
64,rv.moment_generating_func(5f64)?);
65let jp: Vec<Vec<f64>> = (0..10).
66map(|_|{
67(0..10).map(|_| 1f64/200f64).collect()
68}).collect();
69let rvv = RandomVector:: new(vec![rv.clone(),
70rv.clone()],JointProbability::new(jp))?;
71println!("Random Vector Covariance === {:?}"
72,rvv.covariance()?);
73println!("Random Vector Correlation=== {:?}"
74,rvv.correlation()?);
75rvv.get_jointprobability()?.validate_probability();
76println!("\n\n\n\n");
77Ok(())
78
79}
80
81
82impl JointProbability<f64>{
83pub fn new(px: Vec<Vec<f64>>)->Self{
92Self{matrix:px}
93}
94
95pub fn validate_probability(&self){
97println!("Sum of all elements is == {:?}",
98self.matrix.iter().
99fold(0f64,|acc,px2|{
100acc+px2.iter().fold(0f64,|acc2,a|{
101acc2+a
102})
103}));
104}
105}
106
107
108impl RandomVariable{
109 pub fn new(x:Vec<f64>,px:Vec<f64>)->Self{
111Self{variable:x,probability:px}
112}
113
114pub fn mean(&self)->CDHResult<f64>{
125self.expectation_from_func(|a|{
126a
127})
128}
129
130pub fn variance(&self)->CDHResult<f64>{
134let summat:f64 = self.variable.iter().
135zip(self.probability.iter()).
136fold(0.0,|acc,b|{
137acc+b.0.powf(2.0)*b.1
138});
139Ok(
140summat-self.mean()?.powf(2.0)
141)
142}
143
144 pub fn expectation_from_func<H>(&self,h:H)
148->CDHResult<f64>
149where H: Fn(f64)->f64
150{Ok(
151self.variable.iter().zip(self.probability.iter()).
152fold(0.0,|mut a,b|{
153a+=h(*b.0)*b.1;
154a
155}))
156}
157
158pub fn expectation_from_set(&self,hx: &[f64])
160->CDHResult<f64>{
161Ok(self.probability.iter().zip(hx.iter()).
162fold(0.0,|mut a, b|{
163a+= b.0*b.1;
164a
165}))
166}
167
168 pub fn std_deviation(&self)->CDHResult<f64>{
170Ok(self.variance()?.powf(0.5))
171}
172
173pub fn moment_generating_func(&self,t:f64)
177->CDHResult<f64>{
178self.expectation_from_func(|a|{
179let e: f64 = std::f64::consts::E;
180e.powf(t*a)
181})
182}
183
184 #[deprecated(
189 since = "0.1.0",
190 note = "Returns Err. f64 cannot represent complex numbers required \
191 for E[e^(itX)]. Will be implemented when complex number \
192 support is added."
193)]
194pub fn characteristic_func(&self,t:f64)
195->CDHResult<f64>{
196Err("characteristic_func requires complex number support \
202 which is not yet implemented. f64 cannot represent \
203 imaginary numbers.".to_string())
204}
205
206pub fn variance_operator<H>(&self,h:H)
208->CDHResult<f64>
209where
210H:Clone+ Fn(f64)->f64
211{
212self.expectation_from_func(|a|{
213(h(a) - self.expectation_from_func(h.clone())
214.unwrap())
215.powf(2.0)
216})
217
218}
219
220}
221
222
223
224impl RandomVector<f64>{
225pub fn new(x:Vec<RandomVariable>,
230px:JointProbability<f64>)->CDHResult<Self>{
231
232match x{
233a if a.len()<2 => Err("Vec.len() should be >= 2 in RandomCector::new".to_string()),
234a => {
235let len = a[0].variable.len();
236for i in &a[1..]{
237if i.variable.len() !=len{
238return Err("Random Variables should have same size".to_string());
239}
240}
241Ok(Self{variables:a,probability_matrix:px})
242}
243}
244}
245
246 pub fn get_jointprobability(&self)
248->CDHResult<JointProbability<f64>>{
249let jp : Vec<Vec<f64>> = self.variables[1]
250.probability.iter()
251.map(|&px2|{
252self.variables[0].probability.iter().
253map(|&px1| px2*px1).collect()
254}).collect();
255Ok(JointProbability::new(jp))
257}
258
259 pub fn marginal_x2(&self,index:usize)
261->CDHResult<f64>{
262Ok(self.probability_matrix.matrix
263.iter().map(|a| a[index]).sum::<f64>())
264}
265
266 pub fn marginal_x1(&self,index:usize)
268->CDHResult<f64>{
269Ok(self.probability_matrix.matrix[index].iter().sum::<f64>())
270}
271
272 pub fn conditional_x1(&self,
274index_x1:usize,
275index_x2:usize)->CDHResult<f64>{
276
277Ok(
278self.probability_matrix
279.matrix[index_x2][index_x1]/self
280.variables[1].probability[index_x2]
281)
282}
283
284pub fn conditional_x2(&self,
286index_x1:usize,
287index_x2:usize)->CDHResult<f64>{
288Ok(
289self.probability_matrix.
290matrix[index_x2][index_x1]/self.
291variables[0].probability[index_x1]
292)}
293
294pub fn conditional_expectation_x1<H>(&self,
296 index_x2: usize, h: H) -> CDHResult<f64>
297 where
298 H: Fn(f64) -> f64,
299 {
300 let x1_outcomes = &self.variables[0]
301.variable;
302let mut expected_value = 0.0;
304
305 for index_x1 in 0..x1_outcomes.len() {
307 let p_cond = self.conditional_x1(
308index_x1, index_x2)?;
309 if p_cond > 0.0 {
310 expected_value += h(
311x1_outcomes[index_x1]) * p_cond;
312 }
313 }
314
315 Ok(expected_value)
316 }
317
318
319 pub fn conditional_expectation_x2<H>(&self,
321index_x1: usize, h: H) -> CDHResult<f64>
322 where
323 H: Fn(f64) -> f64,
324 {
325 let x2_outcomes = &self.variables[1]
326.variable;
327let mut expected_value = 0.0;
329
330 for index_x2 in 0..x2_outcomes.len() {
332 let p_cond = self.conditional_x2(
333index_x1, index_x2)?;
334 if p_cond > 0.0 {
335 expected_value += h(
336x2_outcomes[index_x2]) * p_cond;
337 }
338 }
339
340 Ok(expected_value)
341 }
342
343
344 pub fn joint_expectation<F>(&self, g: F)
346 -> CDHResult<f64>
347 where
348 F: Fn(f64, f64) -> f64,
349 {
350 let x = &self.variables[0];
351 let y = &self.variables[1];
352 let joint_matrix = &self.probability_matrix.
353matrix;
354
355 let mut expected_value = 0.0;
356
357 for i in 0..x.variable.len() {
358 for j in 0..y.variable.len() {
359 let p_joint = joint_matrix[i][j];
360
361 if p_joint > 0.0 {
362 expected_value += g(
364x.variable[i], y.variable[j]) * p_joint;
365 }
366 }
367 }
368
369 Ok( expected_value)
370 }
371
372
373pub fn covariance(&self) -> CDHResult<f64> {
380 let e_x1 = self.joint_expectation(
383|x, _y| x)?;
384 let e_x2 = self.joint_expectation(
385|_x, y| y)?;
386
387 Ok( self.joint_expectation(
390|x, y| (x - e_x1) * (y - e_x2))?)
391 }
392
393pub fn correlation(&self) -> CDHResult<f64> {
395 let cov = self.covariance()?;
396
397 let e_x1 = self.joint_expectation(
400|x, _y| x)?;
401 let e_x2 = self.joint_expectation(
402|_x, y| y)?;
403
404 let var_x1 = self.joint_expectation(
407|x, _y| (x - e_x1).powi(2))?;
408 let var_x2 = self.joint_expectation(
409|_x, y| (y - e_x2).powi(2))?;
410
411 Ok( if var_x1 > 0.0 && var_x2 > 0.0 {
415 cov / (var_x1.sqrt() * var_x2.sqrt())
416 } else {
417 0.0
418 })
419 }
420}
421
422
423
424impl Probability{
425 pub fn get_probability(sample_space_len: f64,
427favourable_len: f64)->CDHResult<f64>{
428Ok(
429favourable_len/sample_space_len)
430}
431
432pub fn factorial(n:usize)->CDHResult<f64>{
437match n{
438a if a<1=> Ok(1f64),
439a if a >101=>Err("Can't Calculate n greater than 101. n should be between 0 to 101".to_string()),
440_=>{
441Ok((1..=n).fold(1f64,|a,b|{
442a*(b as f64)
443}))
444}
445}
446
447}
448
449
450 pub fn n_p_r(n:usize,r:usize)->CDHResult<f64>{
452match (n,r){
453(0,_) => Err("n cannot be 0".to_string()),
454(_a,0)=>Ok(1f64),
455(a,b) if a==b=>Self::factorial(a),
456(_,_)=>{
457Ok((0..r).fold(1f64,|a,b|{
458a*(n as f64-b as f64)
459}))
460}
461}
462
463
464}
465
466 pub fn n_c_r(n:usize,r:usize)->CDHResult<f64>{
468Ok(Self::n_p_r(n,r)?/Self::factorial(r)?)
469}
470
471 pub fn random_variable<F>(set: &[f64],
473mut f:F)-> CDHResult<Vec<f64>>
474where
475F: FnMut(f64)->f64,
476{Ok(
478set.into_iter().map(|elements|{
479f(*elements)
480}).collect::<Vec<f64>>())
481}
482
483pub fn distribution_function(
485probability_set:&[f64])->CDHResult<Vec<f64>>{
486let mut distri: Vec<f64>= Vec::with_capacity(
487probability_set.len());
488let mut sum = 0f64;
489for &prob in probability_set{
490sum+= prob;
491distri.push(sum);
492}
493Ok(
494distri)
495}
496
497}
498
499
500pub struct Binomial {
502 pub n: usize,
504 pub x: Vec<usize>,
506 pub p: f64,
508}
509
510
511impl Binomial{
512pub fn new(n:usize,x:Vec<usize>,p:f64)->Self{
513Self{
514n,x,p
515}
516}
517
518pub fn get_probability_set(&self)
520->CDHResult<Vec<f64>>{
521Ok(self.x.iter().map(|&i|{
522let ncr = Probability::n_c_r(self.n,i).unwrap();
523let ppx = self.p.powf(i as f64);
524let qpnx = (1f64-self.p).powf((self.n-i) as f64);
525ncr*ppx*qpnx
526}).collect::<Vec<f64>>()
527)}
528
529
530pub fn mean(&self)->CDHResult<f64>{
531Ok(self.n as f64 * self.p)
532}
533pub fn variance(&self)->CDHResult<f64>{
534Ok(self.mean()?*(1f64-self.p))
535}
536}
537
538
539pub struct Geometric {
541 pub x: Vec<usize>,
543 pub p: f64,
545}
546
547
548impl Geometric{
549pub fn new(x:Vec<usize>, p:f64)->Self{
550Self{
551x,p}
552}
553
554pub fn get_probability_set(&self)
556->CDHResult<Vec<f64>>{
557Ok(self.x.iter().map(|&i|{
558let qpnx = (1f64-self.p).powf(i as f64 -1f64);
559self.p*qpnx
560}).collect::<Vec<f64>>()
561)}
562pub fn mean(&self)->CDHResult<f64>{
563Ok(1f64/ self.p)
564}
565pub fn variance(&self)->CDHResult<f64>{
566Ok((1f64-self.p)/(self.p.powf(2f64)))
567}
568}
569
570pub struct Pascal {
572 pub r: usize,
574 pub x: Vec<usize>,
576 pub p: f64,
578}
579
580
581impl Pascal{
582pub fn new(r:usize,x:Vec<usize>,p:f64)->Self{
583Self{r,x,p}
584}
585
586pub fn get_probability_set(&self)
588->CDHResult<Vec<f64>>{
589Ok(self.x.iter().map(|&i|{
590let ncr = Probability::n_c_r(i-1,self.r-1).unwrap();
591let ppx = self.p.powf(i as f64);
592let qpnx = (1f64-self.p).powf((i-self.r) as f64);
593ncr*ppx*qpnx
594}).collect::<Vec<f64>>()
595)}
596
597pub fn mean(&self)->CDHResult<f64>{
598Ok(self.r as f64/self.p)
599}
600
601pub fn variance(&self)->CDHResult<f64>{
602Ok(self.r as f64/(self.p.powf(2f64)))
603}
604}
605
606pub struct HyperGeometric {
608 pub big_n: usize,
610 pub big_t: usize,
612 pub n: usize,
614 pub x: usize,
616}
617
618
619impl HyperGeometric{
620pub fn new(big_n:usize,big_t: usize,
621n:usize,x:usize)->Self{
622Self{big_n,big_t,n,x}
623}
624
625 pub fn get_probability(&self)->CDHResult<f64>{
627
628let big_t_x = Probability::n_c_r(self.big_t,
629self.x)?;
630let big_n_t = Probability::n_c_r(
631self.big_n-self.big_t,
632self.n-self.x)?;
633let big_n_n = Probability::n_c_r(self.big_n,
634self.n)?;
635Ok(big_t_x*big_n_t/big_n_n)
636}
637
638pub fn mean(&self)->CDHResult<f64>{
639Ok(self.n as f64 * self.big_t as f64/self.big_n as f64)
640}
641
642pub fn variance(&self)->CDHResult<f64>{
643let dd = self.big_t as f64;
644let nn= self.big_n as f64;
645let n= self.n as f64;
646Ok(self.mean()?*(1f64 - (dd/nn))*((nn-n)/(nn-1f64)))
647}
648}