equal
deleted
inserted
replaced
505 { |
505 { |
506 if(!m_generator) |
506 if(!m_generator) |
507 { |
507 { |
508 m_generator = new RngStream(); |
508 m_generator = new RngStream(); |
509 } |
509 } |
510 double r = -m_mean*log(m_generator->RandU01()); |
510 while(1) |
511 if (m_bound != 0 && r > m_bound) |
511 { |
512 { |
512 double r = -m_mean*log(m_generator->RandU01()); |
513 return m_bound; |
513 if (m_bound == 0 || r <= m_bound) return r; |
514 } |
514 //otherwise, try again |
515 return r; |
515 } |
516 } |
516 } |
517 |
517 |
518 RandomVariableBase* ExponentialVariableImpl::Copy() const |
518 RandomVariableBase* ExponentialVariableImpl::Copy() const |
519 { |
519 { |
520 return new ExponentialVariableImpl(*this); |
520 return new ExponentialVariableImpl(*this); |
605 if(!m_generator) |
605 if(!m_generator) |
606 { |
606 { |
607 m_generator = new RngStream(); |
607 m_generator = new RngStream(); |
608 } |
608 } |
609 double scale = m_mean * ( m_shape - 1.0) / m_shape; |
609 double scale = m_mean * ( m_shape - 1.0) / m_shape; |
610 double r = (scale * ( 1.0 / pow(m_generator->RandU01(), 1.0 / m_shape))); |
610 while(1) |
611 if (m_bound != 0 && r > m_bound) |
611 { |
612 { |
612 double r = (scale * ( 1.0 / pow(m_generator->RandU01(), 1.0 / m_shape))); |
613 return m_bound; |
613 if (m_bound == 0 || r <= m_bound) return r; |
614 } |
614 //otherwise, try again |
615 return r; |
615 } |
616 } |
616 } |
617 |
617 |
618 RandomVariableBase* ParetoVariableImpl::Copy() const |
618 RandomVariableBase* ParetoVariableImpl::Copy() const |
619 { |
619 { |
620 return new ParetoVariableImpl(*this); |
620 return new ParetoVariableImpl(*this); |
704 if(!m_generator) |
704 if(!m_generator) |
705 { |
705 { |
706 m_generator = new RngStream(); |
706 m_generator = new RngStream(); |
707 } |
707 } |
708 double exponent = 1.0 / m_alpha; |
708 double exponent = 1.0 / m_alpha; |
709 double r = m_mean * pow( -log(m_generator->RandU01()), exponent); |
709 while(1) |
710 if (m_bound != 0 && r > m_bound) |
710 { |
711 { |
711 double r = m_mean * pow( -log(m_generator->RandU01()), exponent); |
712 return m_bound; |
712 if (m_bound == 0 || r <= m_bound) return r; |
713 } |
713 //otherwise, try again |
714 return r; |
714 } |
715 } |
715 } |
716 |
716 |
717 RandomVariableBase* WeibullVariableImpl::Copy() const |
717 RandomVariableBase* WeibullVariableImpl::Copy() const |
718 { |
718 { |
719 return new WeibullVariableImpl(*this); |
719 return new WeibullVariableImpl(*this); |
748 |
748 |
749 /** |
749 /** |
750 * \brief Construct a normal random variable with specified mean and variance |
750 * \brief Construct a normal random variable with specified mean and variance |
751 * \param m Mean value |
751 * \param m Mean value |
752 * \param v Variance |
752 * \param v Variance |
753 * \param b Bound. The NormalVariableImpl is bounded within +-bound. |
753 * \param b Bound. The NormalVariableImpl is bounded within +-bound of the mean. |
754 */ |
754 */ |
755 NormalVariableImpl(double m, double v, double b = INFINITE_VALUE); |
755 NormalVariableImpl(double m, double v, double b = INFINITE_VALUE); |
756 |
756 |
757 NormalVariableImpl(const NormalVariableImpl& c); |
757 NormalVariableImpl(const NormalVariableImpl& c); |
758 |
758 |
763 virtual RandomVariableBase* Copy(void) const; |
763 virtual RandomVariableBase* Copy(void) const; |
764 |
764 |
765 private: |
765 private: |
766 double m_mean; // Mean value of RV |
766 double m_mean; // Mean value of RV |
767 double m_variance; // Mean value of RV |
767 double m_variance; // Mean value of RV |
768 double m_bound; // Bound on value (absolute value) |
768 double m_bound; // Bound on value's difference from the mean (absolute value) |
769 bool m_nextValid; // True if next valid |
769 bool m_nextValid; // True if next valid |
770 double m_next; // The algorithm produces two values at a time |
770 double m_next; // The algorithm produces two values at a time |
771 static bool m_static_nextValid; |
771 static bool m_static_nextValid; |
772 static double m_static_next; |
772 static double m_static_next; |
773 }; |
773 }; |
808 double w = v1 * v1 + v2 * v2; |
808 double w = v1 * v1 + v2 * v2; |
809 if (w <= 1.0) |
809 if (w <= 1.0) |
810 { // Got good pair |
810 { // Got good pair |
811 double y = sqrt((-2 * log(w))/w); |
811 double y = sqrt((-2 * log(w))/w); |
812 m_next = m_mean + v2 * y * sqrt(m_variance); |
812 m_next = m_mean + v2 * y * sqrt(m_variance); |
813 if (fabs(m_next) > m_bound) |
813 //if next is in bounds, it is valid |
|
814 m_nextValid = fabs(m_next-m_mean) <= m_bound; |
|
815 double x1 = m_mean + v1 * y * sqrt(m_variance); |
|
816 //if x1 is in bounds, return it |
|
817 if (fabs(x1-m_mean) <= m_bound) |
814 { |
818 { |
815 m_next = m_bound * (m_next)/fabs(m_next); |
819 return x1; |
816 } |
820 } |
817 m_nextValid = true; |
821 //otherwise try and return m_next if it is valid |
818 double x1 = m_mean + v1 * y * sqrt(m_variance); |
822 else if (m_nextValid) |
819 if (fabs(x1) > m_bound) |
823 { |
820 { |
824 m_nextValid = false; |
821 x1 = m_bound * (x1)/fabs(x1); |
825 return m_next; |
822 } |
826 } |
823 return x1; |
827 //otherwise, just run this loop again |
824 } |
828 } |
825 } |
829 } |
826 } |
830 } |
827 |
831 |
828 RandomVariableBase* NormalVariableImpl::Copy() const |
832 RandomVariableBase* NormalVariableImpl::Copy() const |