diff --git a/include/pops/exponential_power_kernel.hpp b/include/pops/exponential_power_kernel.hpp
index 1abecccf..09473ab6 100644
--- a/include/pops/exponential_power_kernel.hpp
+++ b/include/pops/exponential_power_kernel.hpp
@@ -69,7 +69,7 @@ class ExponentialPowerKernel
     double pdf(double x)
     {
         return (beta / (2 * alpha * std::tgamma(1.0 / beta)))
-               * pow(exp(-x / alpha), beta);
+               * exp(-pow(x / alpha, beta));
     }
 
     /*!
@@ -83,9 +83,13 @@ class ExponentialPowerKernel
         if (x <= 0 || x >= 1) {
             throw std::invalid_argument("icdf: x must be between 0.0 and 1.0");
         }
+        float sign = ((x - 0.5) > 0) ? 1.0 : ((x - 0.5) < 0 ? -1.0 : 0);
+        if (sign == 0) {
+            return 0;
+        }
         GammaKernel gamma_distribution(1.0 / beta, 1.0 / pow(alpha, beta));
         double gamma = gamma_distribution.icdf(2 * std::abs(x - 0.5));
-        return (x - 0.5) * pow(gamma, 1.0 / beta);
+        return sign * pow(gamma, 1.0 / beta);
     }
 };
 
diff --git a/include/pops/gamma_kernel.hpp b/include/pops/gamma_kernel.hpp
index fd3e1073..af57e1c2 100644
--- a/include/pops/gamma_kernel.hpp
+++ b/include/pops/gamma_kernel.hpp
@@ -85,9 +85,9 @@ class GammaKernel
         double beta = 1.0 / theta;
         for (int i = 0; i < alpha; i++) {
             // tgamma = (i-1)! used since c++ has no factorial in std lib
-            sum += pow(beta * x, i) / std::tgamma(i + 1);
+            sum += (1.0 / std::tgamma(i + 1)) * exp(-beta * x) * pow(beta * x, i);
         }
-        return 1 - sum * exp(-beta * x);
+        return 1 - sum;
     }
 
     /*!
@@ -110,7 +110,7 @@ class GammaKernel
         LogNormalKernel lognormal(1);
         double guess = lognormal.icdf(x);
         double check = cdf(guess);
-        double numiterations = 500;  // will need to adjust this
+        double numiterations = 1000;  // will need to adjust this
         double precision = 0.001;  // will need to adjust this
         for (int i = 0; i < numiterations; i++) {
             if (check < (x - precision) || check > (x + precision)) {
@@ -118,16 +118,19 @@ class GammaKernel
                 // if dif is positive guess is greater than needed
                 // if dif is negative guess is less than needed
                 double past_guess = guess;
-                double derivative = (check - x) / pdf(guess);
+                double derivative = dif / pdf(guess);
                 // limit size of next guess
                 guess = std::max(guess / 10, std::min(guess * 10, guess - derivative));
                 check = cdf(guess);
                 // Check if we went to far and need to backtrack
-                for (int j = 0; j < 10; j++) {
-                    if (std::abs(dif) < std::abs(check - x)) {
-                        past_guess = guess;
-                        guess = (guess + past_guess) / 2.0;
-                        check = cdf(guess);
+                int count = 0;
+                bool run = true;
+                while ((std::abs(dif) < std::abs(check - x)) && run) {
+                    guess = (guess + past_guess) / 2.0;
+                    check = cdf(guess);
+                    count++;
+                    if (count > 20) {
+                        run = false;
                     }
                 }
             }
diff --git a/include/pops/lognormal_kernel.hpp b/include/pops/lognormal_kernel.hpp
index e04550db..60abd48e 100644
--- a/include/pops/lognormal_kernel.hpp
+++ b/include/pops/lognormal_kernel.hpp
@@ -72,7 +72,7 @@ class LogNormalKernel
         if (x == 0) {
             return 0;
         }
-        return (1 / x) * (1 / (sigma * sqrt(2 * M_PI)))
+        return (1 / (x * sigma * sqrt(2 * M_PI)))
                * exp(-(pow(log(x), 2)) / (2 * pow(sigma, 2)));
     }
 
@@ -93,9 +93,14 @@ class LogNormalKernel
         double y = (2 * x) - 1;
         float sign = (y < 0) ? -1.0f : 1.0f;
         //  0.147 used for a relative error of about 2*10^-3
-        float b = 2 / (M_PI * 0.147) + 0.5f * log(1 - pow(y, 2));
+        // float b = 2 / (M_PI * 0.147) + 0.5f * log(1 - pow(y, 2));
+        // double inverf =
+        //    (sign * sqrt(-b + sqrt(pow(b, 2) - (1 / (0.147) * log(1 - pow(y, 2))))));
+        double a = 0.140012;
+        double t = 2.0 / (M_PI * a);
+        double l = log(1 - pow(y, 2));
         double inverf =
-            (sign * sqrt(-b + sqrt(pow(b, 2) - (1 / (0.147) * log(1 - pow(y, 2))))));
+            sign * sqrt(sqrt(pow(t + (l / 2.0), 2) - (l / a)) - (t + (l / 2.0)));
         return exp(sqrt(2 * pow(sigma, 2)) * inverf);
     }
 };
diff --git a/include/pops/normal_kernel.hpp b/include/pops/normal_kernel.hpp
index aa59d913..e4c8bb60 100644
--- a/include/pops/normal_kernel.hpp
+++ b/include/pops/normal_kernel.hpp
@@ -87,11 +87,11 @@ class NormalKernel
         //  approximation for inverse error function
         double y = (2 * x) - 1;
         float sign = (y < 0) ? -1.0f : 1.0f;
-        //  0.147 used for a relative error of about 2*10^-3
-        float b = 2.0 / (M_PI * 0.147) + 0.5f * log(1 - pow(y, 2));
+        double a = 0.140012;
+        double t = 2.0 / (M_PI * a);
+        double l = log(1 - pow(y, 2));
         double inverf =
-            (sign * sqrt(-b + sqrt(pow(b, 2) - (1.0 / (0.147) * log(1 - pow(y, 2))))));
-
+            sign * sqrt(sqrt(pow(t + (l / 2.0), 2) - (l / a)) - (t + (l / 2.0)));
         return sigma * std::sqrt(2) * inverf;
     }
 };
diff --git a/include/pops/power_law_kernel.hpp b/include/pops/power_law_kernel.hpp
index 138feb58..35f7eb78 100644
--- a/include/pops/power_law_kernel.hpp
+++ b/include/pops/power_law_kernel.hpp
@@ -31,12 +31,12 @@ using std::pow;
 class PowerLawKernel
 {
 protected:
-    double xmin;
     double alpha;
+    double xmin;
     std::uniform_real_distribution<double> distribution;
 
 public:
-    PowerLawKernel(double xm, double a) : xmin(xm), alpha(a), distribution(0.0, 1.0)
+    PowerLawKernel(double a, double xm) : alpha(a), xmin(xm), distribution(0.0, 1.0)
     {
         //         if (alpha < 1.0) {
         //             throw std::invalid_argument("alpha must be greater than or equal
@@ -92,7 +92,7 @@ class PowerLawKernel
         if (x <= 0 || x >= 1) {
             throw std::invalid_argument("icdf: x must be between 0 and 1.0");
         }
-        return xmin * pow(x, (1.0 / (-alpha + 1.0)));
+        return pow(x / xmin, (-alpha + 1.0));
     }
 };
 
diff --git a/include/pops/radial_kernel.hpp b/include/pops/radial_kernel.hpp
index 9e61b35f..82341878 100644
--- a/include/pops/radial_kernel.hpp
+++ b/include/pops/radial_kernel.hpp
@@ -184,7 +184,7 @@ class RadialDispersalKernel
         else if (dispersal_kernel_type_ == DispersalKernelType::Gamma) {
             distance = std::abs(gamma_distribution.random(generator));
         }
-        else if (dispersal_kernel_type_ == DispersalKernelType::Logistic) {
+        else if (dispersal_kernel_type_ == DispersalKernelType::ExponentialPower) {
             distance = std::abs(exponential_power_distribution.random(generator));
         }
         else if (dispersal_kernel_type_ == DispersalKernelType::Logistic) {
diff --git a/tests/test_deterministic.cpp b/tests/test_deterministic.cpp
index 9b5fd6c8..cd6f1833 100644
--- a/tests/test_deterministic.cpp
+++ b/tests/test_deterministic.cpp
@@ -26,6 +26,9 @@
 #include <pops/deterministic_kernel.hpp>
 #include <pops/cauchy_kernel.hpp>
 #include <pops/exponential_kernel.hpp>
+#include <pops/exponential_power_kernel.hpp>
+#include <pops/lognormal_kernel.hpp>
+#include <pops/gamma_kernel.hpp>
 #include <pops/simulation.hpp>
 
 using std::string;
@@ -57,6 +60,7 @@ int test_with_cauchy_deterministic_kernel()
     double reproductive_rate = 2;
     bool generate_stochasticity = false;
     bool establishment_stochasticity = false;
+    bool movement_stochasticity = false;
     // We want everything to establish.
     double establishment_probability = 1;
     // Cauchy
@@ -67,7 +71,8 @@ int test_with_cauchy_deterministic_kernel()
         model_type_from_string("SI"),
         0,
         generate_stochasticity,
-        establishment_stochasticity);
+        establishment_stochasticity,
+        movement_stochasticity);
     simulation.generate(
         dispersers,
         infected,
@@ -141,6 +146,7 @@ int test_with_exponential_deterministic_kernel()
     double reproductive_rate = 2;
     bool generate_stochasticity = false;
     bool establishment_stochasticity = false;
+    bool movement_stochasticity = false;
     // We want everything to establish.
     double establishment_probability = 1;
     // Exponential
@@ -151,7 +157,8 @@ int test_with_exponential_deterministic_kernel()
         model_type_from_string("SI"),
         0,
         generate_stochasticity,
-        establishment_stochasticity);
+        establishment_stochasticity,
+        movement_stochasticity);
     s2.generate(
         dispersers,
         infected,
@@ -226,6 +233,7 @@ int test_with_weibull_deterministic_kernel()
     double reproductive_rate = 2;
     bool generate_stochasticity = false;
     bool establishment_stochasticity = false;
+    bool movement_stochasticity = false;
     // We want everything to establish.
     double establishment_probability = 1;
     Simulation<Raster<int>, Raster<double>> s2(
@@ -235,7 +243,8 @@ int test_with_weibull_deterministic_kernel()
         model_type_from_string("SI"),
         0,
         generate_stochasticity,
-        establishment_stochasticity);
+        establishment_stochasticity,
+        movement_stochasticity);
     s2.generate(
         dispersers,
         infected,
@@ -309,6 +318,7 @@ int test_with_log_normal_deterministic_kernel()
     double reproductive_rate = 2;
     bool generate_stochasticity = false;
     bool establishment_stochasticity = false;
+    bool movement_stochasticity = false;
     // We want everything to establish.
     double establishment_probability = 1;
     Simulation<Raster<int>, Raster<double>> s2(
@@ -318,7 +328,8 @@ int test_with_log_normal_deterministic_kernel()
         model_type_from_string("SI"),
         0,
         generate_stochasticity,
-        establishment_stochasticity);
+        establishment_stochasticity,
+        movement_stochasticity);
     s2.generate(
         dispersers,
         infected,
@@ -393,6 +404,7 @@ int test_with_normal_deterministic_kernel()
     double reproductive_rate = 2;
     bool generate_stochasticity = false;
     bool establishment_stochasticity = false;
+    bool movement_stochasticity = false;
     // We want everything to establish.
     double establishment_probability = 1;
     Simulation<Raster<int>, Raster<double>> s2(
@@ -402,7 +414,8 @@ int test_with_normal_deterministic_kernel()
         model_type_from_string("SI"),
         0,
         generate_stochasticity,
-        establishment_stochasticity);
+        establishment_stochasticity,
+        movement_stochasticity);
     s2.generate(
         dispersers,
         infected,
@@ -476,6 +489,7 @@ int test_with_hyperbolic_secant_deterministic_kernel()
     double reproductive_rate = 2;
     bool generate_stochasticity = false;
     bool establishment_stochasticity = false;
+    bool movement_stochasticity = false;
     // We want everything to establish.
     double establishment_probability = 1;
     Simulation<Raster<int>, Raster<double>> s2(
@@ -485,7 +499,8 @@ int test_with_hyperbolic_secant_deterministic_kernel()
         model_type_from_string("SI"),
         0,
         generate_stochasticity,
-        establishment_stochasticity);
+        establishment_stochasticity,
+        movement_stochasticity);
     s2.generate(
         dispersers,
         infected,
@@ -548,8 +563,8 @@ int test_with_power_law_deterministic_kernel()
     Raster<double> weather_coefficient = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
     std::vector<unsigned> movement_schedule = {1, 1};
 
-    Raster<int> expected_mortality_tracker = {{10, 0, 0}, {0, 10, 0}, {0, 0, 2}};
-    Raster<int> expected_infected = {{15, 0, 0}, {0, 15, 0}, {0, 0, 4}};
+    Raster<int> expected_mortality_tracker = {{9, 1, 0}, {0, 9, 0}, {0, 0, 2}};
+    Raster<int> expected_infected = {{14, 1, 0}, {0, 14, 0}, {0, 0, 4}};
 
     Raster<int> dispersers(infected.rows(), infected.cols());
     std::vector<std::tuple<int, int>> outside_dispersers;
@@ -561,6 +576,7 @@ int test_with_power_law_deterministic_kernel()
     double reproductive_rate = 2;
     bool generate_stochasticity = false;
     bool establishment_stochasticity = false;
+    bool movement_stochasticity = false;
     // We want everything to establish.
     double establishment_probability = 1;
     Simulation<Raster<int>, Raster<double>> s2(
@@ -570,7 +586,8 @@ int test_with_power_law_deterministic_kernel()
         model_type_from_string("SI"),
         0,
         generate_stochasticity,
-        establishment_stochasticity);
+        establishment_stochasticity,
+        movement_stochasticity);
     s2.generate(
         dispersers,
         infected,
@@ -600,7 +617,7 @@ int test_with_power_law_deterministic_kernel()
         deterministicKernel,
         suitable_cell,
         establishment_probability);
-    if (outside_dispersers.size() != 0) {
+    if (outside_dispersers.size() != 1) {
         cout << "Deterministic Kernel PowerLaw: There are outside_dispersers ("
              << outside_dispersers.size() << ") but there should be 0\n";
         return 1;
@@ -643,6 +660,7 @@ int test_with_logistic_deterministic_kernel()
     double reproductive_rate = 2;
     bool generate_stochasticity = false;
     bool establishment_stochasticity = false;
+    bool movement_stochasticity = false;
     // We want everything to establish.
     double establishment_probability = 1;
     Simulation<Raster<int>, Raster<double>> s2(
@@ -652,7 +670,8 @@ int test_with_logistic_deterministic_kernel()
         model_type_from_string("SI"),
         0,
         generate_stochasticity,
-        establishment_stochasticity);
+        establishment_stochasticity,
+        movement_stochasticity);
     s2.generate(
         dispersers,
         infected,
@@ -725,6 +744,7 @@ int test_with_gamma_deterministic_kernel()
     double reproductive_rate = 2;
     bool generate_stochasticity = false;
     bool establishment_stochasticity = false;
+    bool movement_stochasticity = false;
     // We want everything to establish.
     double establishment_probability = 1;
     Simulation<Raster<int>, Raster<double>> s2(
@@ -734,7 +754,8 @@ int test_with_gamma_deterministic_kernel()
         model_type_from_string("SI"),
         0,
         generate_stochasticity,
-        establishment_stochasticity);
+        establishment_stochasticity,
+        movement_stochasticity);
     s2.generate(
         dispersers,
         infected,
@@ -808,6 +829,7 @@ int test_with_exponential_power_deterministic_kernel()
     double reproductive_rate = 2;
     bool generate_stochasticity = false;
     bool establishment_stochasticity = false;
+    bool movement_stochasticity = false;
     // We want everything to establish.
     double establishment_probability = 1;
     Simulation<Raster<int>, Raster<double>> s2(
@@ -817,7 +839,8 @@ int test_with_exponential_power_deterministic_kernel()
         model_type_from_string("SI"),
         0,
         generate_stochasticity,
-        establishment_stochasticity);
+        establishment_stochasticity,
+        movement_stochasticity);
     s2.generate(
         dispersers,
         infected,
@@ -943,6 +966,77 @@ int test_exponential_distribution_functions()
     return 0;
 }
 
+int test_gamma_distribution_functions()
+{
+    cout << "Gamma\n";
+    double x = 0.1;
+    for (double a = 0.1; a < 30; a += 5) {
+        for (double t = 0.1; t < 30; t += 5) {
+            x = 0.1;
+            for (double icdf_x = 0.1; icdf_x < 1; icdf_x += 0.1) {
+                // testing gamma pdf & icdf
+                cout << "alpha = " << a << " theta = " << t << " x = " << x << '\n';
+                GammaKernel gamma(a, t);
+                double probability = gamma.pdf(x);
+                cout << " pdf = " << probability << '\n';
+                double icdf = -11;
+                try {
+                    icdf = gamma.icdf(icdf_x);
+                }
+                catch (const std::invalid_argument& ia) {
+                    std::cerr << "Invalid argument: " << ia.what() << "x = " << x
+                              << " alpha = " << a << " theta = " << t << '\n';
+                }
+                cout << " icdf = " << icdf << "\n";
+                x += 5;
+            }
+        }
+    }
+    return 1;
+}
+
+int test_exponential_power_distribution_functions()
+{
+    cout << "Exponential Power\n";
+    double x = 0;
+    for (double a = 1; a < 5; a += 1) {
+        for (double b = 0.5; b < 5; b += 1) {
+            x = 0;
+            for (double icdf_x = 0.1; icdf_x < 1; icdf_x += 0.1) {
+                // testing gamma pdf & icdf
+                ExponentialPowerKernel ep(a, b);
+                cout << "alpha = " << a << " beta = " << b << " x = " << x;
+                double probability = ep.pdf(x);
+                cout << " pdf = " << probability << " icdf x = " << icdf_x;
+                double icdf = ep.icdf(icdf_x);
+                cout << " icdf = " << icdf << "\n";
+                x += 0.50;
+            }
+        }
+    }
+    return 1;
+}
+
+int test_log_normal_distribution_functions()
+{
+    cout << "Log Normal\n";
+    double x = 0;
+    for (double s = 0.1; s < 10; s += 1) {
+        x = 0;
+        for (double icdf_x = 0.1; icdf_x < 1; icdf_x += 0.1) {
+            // testing gamma pdf & icdf
+            LogNormalKernel ln(s);
+            cout << "sigma = " << s << " x = " << x;
+            double probability = ln.pdf(x);
+            cout << " pdf = " << probability << " icdf x = " << icdf_x;
+            double icdf = ln.icdf(icdf_x);
+            cout << " icdf = " << icdf << "\n";
+            x += 0.25;
+        }
+    }
+    return 1;
+}
+
 int main()
 {
     int ret = 0;
@@ -959,6 +1053,9 @@ int main()
     ret += test_with_logistic_deterministic_kernel();
     ret += test_with_gamma_deterministic_kernel();
     ret += test_with_exponential_power_deterministic_kernel();
+    // ret += test_gamma_distribution_functions();
+    // ret += test_exponential_power_distribution_functions();
+    // ret += test_log_normal_distribution_functions();
 
     std::cout << "Test deterministic number of errors: " << ret << std::endl;
     return ret;
diff --git a/tests/test_random.cpp b/tests/test_random.cpp
index ef0f1def..1dd5e9a4 100644
--- a/tests/test_random.cpp
+++ b/tests/test_random.cpp
@@ -61,10 +61,11 @@ int test_power_law_rng()
 
     double alpha = 1.5;
     double xmin = 0.01;
+    PowerLawKernel power_law_distribution(alpha, xmin);
     std::uniform_real_distribution<double> distribution(0.1, 1.0);
     for (int i = 0; i < 10000; i++) {
         double x = distribution(generator);
-        double y = pow(x, (1.0 / (-alpha + 1.0))) * xmin;
+        double y = power_law_distribution.icdf(x);
         file << x << "," << y << ",\n";
     }
 
@@ -81,6 +82,7 @@ int test_hyperbolic_secant_rng()
     // std::cout << "Hyperbolic secant random: " << logsech.random(generator) << "\n";
     // Pulling code out of hyperbolic_secant_kernel to test here to be able to save the
     // value produced by the rng
+    HyperbolicSecantKernel hyperbolic_secant_distribution(sigma_);
     std::ofstream file;
     file.open("testing_secant.csv");
     file << "x,y,\n";
@@ -88,10 +90,9 @@ int test_hyperbolic_secant_rng()
     std::uniform_real_distribution<double> distribution(0.0, 1.0);
     for (int i = 0; i < 10000; i++) {
         double x = distribution(generator);
-
         // get random value from a uniform distribution and use it
         // to get a random value from the distribution
-        double y = ((log(tan((x * M_PI) / 2.0)) * (2.0 * sigma_)) / M_PI);
+        double y = hyperbolic_secant_distribution.icdf(x);
         file << x << "," << y << ",\n";
     }
     file.close();
@@ -101,6 +102,7 @@ int test_logistic_rng()
 {
     std::default_random_engine generator;
     double s_ = 2.0;
+    LogisticKernel logistic_distribution(s_);
     std::ofstream file;
     file.open("testing_logistic.csv");
     file << "x,y,\n";
@@ -111,7 +113,7 @@ int test_logistic_rng()
 
         // get random value from a uniform distribution and use it
         // to get a random value from the distribution
-        double y = s_ * log(x / (1 - x));
+        double y = logistic_distribution.icdf(x);
         file << x << "," << y << ",\n";
     }
     file.close();
@@ -123,6 +125,7 @@ int test_exponential_power_rng()
     std::default_random_engine generator;
     double alpha_ = 1.5;
     double beta_ = 0.5;
+    ExponentialPowerKernel exponential_power_distribution(alpha_, beta_);
     std::ofstream file;
     file.open("testing_exponential_power.csv");
     file << "x,y,\n";
@@ -130,9 +133,26 @@ int test_exponential_power_rng()
     std::uniform_real_distribution<double> distribution(0.0, 1.0);
     for (int i = 0; i < 10000; i++) {
         double x = distribution(generator);
-        GammaKernel gamma_distribution(1.0 / beta_, 1.0 / pow(alpha_, beta_));
-        double gamma = gamma_distribution.icdf(2 * std::abs(x - 0.5));
-        double y = (x - 0.5) * pow(gamma, 1.0 / beta_);
+        double y = exponential_power_distribution.icdf(x);
+        file << x << "," << y << ",\n";
+    }
+    file.close();
+    return 0;
+}
+
+int test_gamma()
+{
+    std::default_random_engine generator;
+    double alpha_ = 1.5;
+    double theta_ = 0.5;
+    GammaKernel gamma_distribution(alpha_, theta_);
+    std::ofstream file;
+    file.open("testing_gamma.csv");
+    file << "x,y,\n";
+    std::uniform_real_distribution<double> distribution(0.0, 1.0);
+    for (int i = 0; i < 10000; i++) {
+        double x = distribution(generator);
+        double y = gamma_distribution.icdf(x);
         file << x << "," << y << ",\n";
     }
     file.close();
@@ -146,6 +166,7 @@ int main()
     ret += test_hyperbolic_secant_rng();
     ret += test_logistic_rng();
     ret += test_exponential_power_rng();
+    ret += test_gamma();
 
     std::cout << "Test deterministic number of errors: " << ret << std::endl;
     return ret;