DEV Community: Rama Reksotinoyo

Handling Happy Numbers: Building a Set-like Data Structure

Rama Reksotinoyo — Sun, 02 Nov 2025 14:05:13 +0000

Today, as I was tackling some LeetCode problems (easy ones, of course :D), I stumbled upon an interesting one called Happy Number.

According to wikipedia, happy number is a set of number theory, a happy number is a number which eventually reaches 1 when the number is replaced by the sum of the square of each digit.

I.e. images below in wikipedia:

Creating set-like data structure

To check whether a value has already been processed (to avoid an infinite loop), I used a Set data structure. Since I’ve worked with Python before, I’m familiar with the built-in set type, which only allows unique values.

But in Go, there’s no built-in Set type, so I had to create one from scratch and implement some methods to mimic the same behavior.

type Set struct{
    elements map[string] struct{}
}

func NewSet() *Set{
    return &Set{
        elements: make(map[string]struct{}),
    }
}

func (s *Set) Add(value string){
    s.elements[value] = struct{}{}
}

func (s *Set) Contains(value string) bool {
    _, found := s.elements[value]
    return found
}

I define a Set struct with a map where:

Key: string - the actual value i want to store
Value: struct{} - an empty struct that uses zero bytes of memory

Why struct{}? Because i only care about the keys (whether a value exists), not the values. Just try to be as efficient as i can.
In addition, I also defined the NewSet method as a constructor and created two simple methods, Add & Contains, so that the Set I defined behaves like a typical Set.

Logic to solve happy number proble

func sumSquares(n int) int {
    // Convert numbers to strings so they can be accessed character by character
    numStr := strconv.Itoa(n)
    sum := 0

    for _, char := range numStr {
        // Convert characters back to numbers
        digit := int(char - '0')
        sum += digit * digit
    }

    return sum
}

func isHappy(n int) bool {
    seen := NewSet()

    for n != 1 {
        key := strconv.Itoa(n)
        if seen.Contains(key) {
            return false 
        }
        seen.Add(key)
        n = sumSquares(n)
    }
    return true
}

The sumSquares function calculates the sum of the squares of each digit of the number n by: converting n to a string, taking each character, converting it to a number, squaring it, and then adding them together.

The isHappy function is the main function in this program. First, it initializes a Set. Then, it loops as long as the value of n is not equal to 1. Inside the loop, the program checks whether the current number (converted to a string) already exists in the Set. If it does, this indicates a cycle has been detected, and the function returns false. If it doesn't exist yet, the number is added to the Set, and n is updated with the result from the sumSquares function. When n reaches 1, the loop terminates and the function returns true, meaning the number is a happy number.

Run the test

func TestIsHappy(t *testing.T) {
    tests := []struct {
        input    int
        expected bool
    }{
        {1, true},
        {19, true},
        {7, true},
        {2, false},
        {4, false},
        {100, true},
        {24323423, false},
    }

    for _, tt := range tests {
        result := isHappy(tt.input)
        if result != tt.expected {
            t.Errorf("isHappy(%d) = %v; want %v", tt.input, result, tt.expected)
        }
    }
}

I write some testcases (using table test), testing several numerical values that fall into the categories of happy numbers and non-happy numbers.

Benchmark test

func Benchmark(b *testing.B) {
    numbers := []int{1, 7, 19, 100, 12345, 999999}
    for i := 0; i < b.N; i++ {
        for _, num := range numbers {
            isHappy(num)
        }
    }
}

My isHappy() function takes an average of 1.547 nanoseconds to execute once on a single number (in a loop on numbers), and Go runs it hundreds of thousands of times to get an accurate average.
I know this isn't the best result possible, but there may be some improvement in the future.

Escape analysis

Let's talk about a subtle performance trade-off in our implementation. When i run Go's escape analysis tool, it tells us something interesting:

❯ go build -gcflags=-m happy_number/main.go
happy_number/main.go:13:9: &Set{...} escapes to heap

That innocent-looking & in our NewSet() function has consequences.

It feels natural—after all, i'im creating a new Set and want to use it throughout our function. But that pointer return type triggers something called "heap escape."

What's Actually Happening?

Every time i call NewSet():

Go allocates memory on the heap (slow), because of returning pointer
The garbage collector later has to clean it up

Before You Go...

Okay, confession time: I know this code could be better. 😅

The sumSquares function doing all those string conversions? Yeah, a bit overkill. And returning a pointer for a tiny Set struct? Probably not necessary.

Way to Avoid the Pitfall of Cherry-Picking Data

Rama Reksotinoyo — Sat, 02 Aug 2025 07:21:12 +0000

I often cherry-pick samples from a population, even though there are many talented statisticians out there who have developed appropriate sampling methods for sampling needs. To avoid this fallacy, I am forced to learn about probability sampling.

In this post, I'll talk about implementing cluster sampling in C. Cluster sampling is one of those handy statistical sampling methods for pulling samples from a population. There are other cool probability sampling techniques too—like simple random sampling, stratified random sampling, and others—but let’s focus on cluster sampling for now.

Why bother with cluster sampling? Well, when you're dealing with a huge population, building a full sampling frame can be a pain—it’s time-consuming, expensive, and just overkill. That’s where probability sampling comes to the rescue!

And if your data is geographically distributed and naturally forms clusters (like neighborhoods, schools, or districts), then cluster sampling is a no-brainer. It’s efficient, practical, and saves you a ton of effort.

So, let’s dive into the game!

Steps to do the cluster sampling:

Group the member of population into 𝑁 clusters.
Do an SRS (simple random samling) to choose 𝑛 from 𝑁 clusters.
Then, estimate the population using all samples in the selected clusters.

Assume this dataset comes from AirAvata, a startup company that provides delivery services. They have expertise in air delivery. They plan to expand their services to a city, let's call it City A. However, AirAvata needs to estimate the appropriate delivery rate for the residents of City A.

Objectives:

Estimate the average people spend per month on goods delivery in city A.
Estimate the 95% CI for the average people spend per month on goods delivery.

Asumtion:

There are 451 districts in city A.

Considering that City A has been divided into districts, I am thinking of conducting cluster sampling. Here is a screenshot of the dataset.

District : the sampled District
People : the number of people in each districts
Total Spend : total expense in each districts (in Rupiah)

Based on the dataset above, we have information below:

Number of all district in city A = 415 districts
Number of sampled district in city A = 20 districts

To estimate the average people spend per month on goods delivery, we use:

I took the above equation from the statistics program material by pacmann

In the code below, i represent the dataset with an object, then I implement the above equation in several functions.

typedef struct 
{
    int district;
    int people;
    double total_spend;
} SurveyData;

double calculate_sc_squared(SurveyData* data, int count, double r) {
    double sum = 0.0;
    for (int i = 0; i < count; i++) {
        double term = data[i].total_spend - (r * data[i].people);
        sum += pow(term, 2);
    }
    return sum / (count - 1);
}

double calculate_variance_tau_hat(int N, int n, double sc_squared) {
    return N * (N - n) * (sc_squared / n);
}

// Modified calculate_and_display_metrics function
void calculate_and_display_metrics(SurveyData* data, int count) {
    int total_people = 0;
    double total_spend = 0.0;
    double average_spend;

    // Calculate totals
    for (int i = 0; i < count; i++) {
        total_people += data[i].people;
        total_spend += data[i].total_spend;
    }

    printf("Total people across all districts: %d\n", total_people);
    printf("Total spend across all districts: %.2f\n", total_spend);

    if (total_people > 0) {
        average_spend = total_spend / (double)total_people;
        printf("\nAverage monthly spend per person on goods delivery: %.2f\n", average_spend);

        // Calculate variance components
        double sc_squared = calculate_sc_squared(data, count, average_spend);
        double var_y_total_est = TOTAL_DISTRICTS * (TOTAL_DISTRICTS - count) * (sc_squared / count);

        printf("\nVariance Calculations:\n");
        printf("Sampled cluster variance (sc²): %.4f\n", sc_squared);
        printf("Variance of total pop. est.: (Rp %.0f)^2\n", sqrt(var_y_total_est));

        // Calculate mean people per district and total population estimate
        double M_bar = (double)total_people / count;
        double M_tot_est = TOTAL_DISTRICTS * M_bar;

        printf("\nPopulation Estimates:\n");
        printf("Mean people per district: %.0f\n", M_bar);
        printf("Estimated total population in city A: %.0f\n", M_tot_est);

        // Variance of mean per person estimator
        double var_y_mean_est = var_y_total_est / pow(M_tot_est, 2);
        printf("Variance of mean per person est.: (Rp %.0f)^2\n", sqrt(var_y_mean_est));

        // Confidence interval calculations
        printf("\nConfidence Interval:\n");
        printf("z-stat (95%% CI): %.2f\n", Z_95_CI);
        double d = Z_95_CI * sqrt(var_y_mean_est);
        printf("Margin of error (d): Rp %.0f\n", d);
        printf("95%% Confidence Interval: Rp %.2f ± Rp %.0f\n", 
               average_spend, d);
        printf("                      (Rp %.0f - Rp %.0f)\n",
               average_spend - d, average_spend + d);
    } else {
        printf("\nError: No people data available for calculation\n");
    }
}

Here's the program that have been compiled to executable files:

District,People,Total Spend

412,171,64050000
324,185,64500000
182,363,123900000
77,148,55200000
418,257,87450000
185,311,107250000
124,244,74400000
320,416,129750000
332,330,104700000
16,305,100500000
352,272,94800000
417,146,49050000
10,195,65700000
32,238,86250000
94,378,126300000
420,239,78900000
214,326,111600000
200,197,69750000
308,182,66750000
268,318,96750000
data: 20 

Number of all districts in city A: 415
Number of sampled districts in city A: 20
Total records: 21

Total people across all districts: 5221
Total spend across all districts: 1757550000.00

Average monthly spend per person on goods delivery: 336630.91

Variance Calculations:
Sampled cluster variance (sc²): 26330932965053.5312
Variance of total pop. est.: (Rp 464558833)^2

Population Estimates:
Mean people per district: 261
Estimated total population in city A: 108336
Variance of mean per person est.: (Rp 4288)^2

Confidence Interval:
z-stat (95% CI): 1.96
Margin of error (d): Rp 8405
95% Confidence Interval: Rp 336630.91 ± Rp 8405
                      (Rp 328226 - Rp 345036)

Full implementation in C, was posted in github Here

Mean of number of people per district: 261
Estimation of total number of people in city A: 108336
The variance of mean ssu. est.: (Rp 4288)^2

With 95% confidence, we estimate that the true average monthly spending on goods delivery in City A falls between Rp 328,226 and Rp 345,036.

Reference:
Pacmann statistic course (RECOMENDED)
GeeksForGeeks

Enhancing Image Quality with Uniform Pixel Intensity Distribution

Rama Reksotinoyo — Sat, 16 Mar 2024 18:53:19 +0000

A week ago i saw my deprecated materials of image processing when i was in 5 semester on my bachelor degree, and was like, wow this method is gorgeous (in my eyes). I opened the code i wrote on that momment, and i didnt found anything cause i code on matlab, whose method is brutally abstracted (with all due respect with matlab coder).

What is histogram equalization
Histogram equalization is one of the technique for increase the contrass of the image, with this the distribution of the data is made uniform or the data is made to have an even uniform distribution. The mathematical calculation is by changing the degree of grayness of a pixel with a new degree of grayness with a transformation function

Goal

Image from mathwork by matlab

Step

First that i gonna do is to find the distribution of each pixels.
Find the probability density function.
Find the cumulative distribution function.
Calculate the cumulative distribution function with the maximum value, in this case the value is 2255 (8-bit).
Find the histogram equalizaiton level.
Decode the data and transform to image.

code
To make my objectives easier i'm using golang, to get there im using Image library as golang standaed library for image handling. If you are curious why I use go, because I am deepening my ability to use golang language. There may be some of my future writing content about rust because of this news.

Take a seat and get your coffee, this might be a bit boring.

type ImageData struct {
    Pixels [][]color.Gray
    Width  int
    Height int
}

type Histogram struct {
    Nk          map[uint8]int
    Pdf         map[uint8]float64
    Cdf         map[uint8]float64
    CdfMultiply map[uint8]float64
    HistEQ      map[uint8]int
}

type ImageProcessor struct {
}

First, I set up a few objects that will represent the image and the method or step used, in this case I am using a grayscale image.

func (ip *ImageProcessor) LoadImage(filepath string) (*ImageData, error) {
    // gray, err := ioutils.LoadImage(filepath)

    file, err := os.Open(filepath)
    if err != nil {
        return nil, fmt.Errorf("failed to load file: %v", err)
    }
    defer file.Close()

    img, err := png.Decode(file)
    if err != nil {
        return nil, fmt.Errorf("failed to open png format: %v", err)
    }

    bounds := img.Bounds()
    width, height := bounds.Max.X, bounds.Max.Y

    pixels := make([][]color.Gray, height)
    for y := 0; y < height; y++ {
        pixels[y] = make([]color.Gray, width)
        for x := 0; x < width; x++ {
            pixels[y][x] = color.GrayModel.Convert(img.At(x, y)).(color.Gray)
        }
    }

    if err != nil {
        return nil, err
    }

    new_width := len(pixels[0])
    new_height := len(pixels)

    return &ImageData{
        Pixels: pixels,
        Width:  new_width,
        Height: new_height,
    }, nil
}

The LoadImage function, part of the ImageProcessor type, serves the purpose of reading an image file. It takes a file path as input, and then retrieve the image data, and calculates the width and height of the image. In case of any errors during the loading process, it returns nil along with the encountered error. Upon successful loading, it constructs and returns an ImageData struct containing the pixel information (in grayscale), as well as the width and height of the loaded image. In essence, this function facilitates the loading of image data from a specified file, offering crucial details about the pixel composition and dimensions of the image.

func (ip *ImageProcessor) Nk(pix []uint8) map[uint8]int {
    count := make(map[uint8]int)
    for i := 0; i <= 255; i++ {
        count[uint8(i)] = 0
    }
    for _, v := range pix {
        count[v]++
    }
    return count
}

Furthermore, the Nk method aims to find the distribution of each pixel in the input image. Yes, it is just a function to find the number of each pixel, nothing special.

Then, the function above is the step to find pdf (probability density function). The pdf is useful to help calculate the probability that the value of x is in an interval. What needs to be underlined, pdf is not a probability, because it only looks at the height of the area under the curve on a distribution curve.

func (ip *ImageProcessor) Pdf(pix map[uint8]int, nPix int) map[uint8]float64 {
    prob := make(map[uint8]float64)

    var maxPixValue int
    for _, v := range pix {
        if v > maxPixValue {
            maxPixValue = v
        }
    }

    for i := 0; i <= maxPixValue; i++ {
        pixel := uint8(i)
        v, exists := pix[pixel]
        if exists {
            prob[pixel] = float64(v) / float64(nPix)
        } else {
            prob[pixel] = 0.0
        }
    }
    return prob
}

Method above is designed to calculate the probability density function (PDF) of pixel intensities in a grayscale image. It takes two parameters: pix, a map representing the count of occurrences for each pixel intensity, and nPix, the total number of pixels in the image. The function iterates through possible pixel intensity values (from 0 to a specified maximum value, in this case, 215), calculating the probability of each intensity's occurrence in the image. If a given intensity is present in the pix map, it computes the probability as the ratio of its count to the total number of pixels (nPix). If the intensity is not found in the map, indicating no occurrences, the probability is set to 0.0. The resulting probabilities are stored in a map and returned, providing a representation of the likelihood distribution of pixel intensities in the grayscale image.

Next is the calculation of the cumulative distribution function (cdf). If pdf is the frequency or height of a value on a curve, then cdf is used to calculate the probability of the integral of the pdf function, in this case cdf in discrete space.

func (ip *ImageProcessor) Cdf(prob map[uint8]float64) map[uint8]float64 {
    sk := make(map[uint8]float64)
    var temp float64

    for i := 0; i <= len(prob)-1; i++ {
        pixel := uint8(i)
        if i == 0 {
            sk[pixel] = prob[pixel]
        } else {
            temp = sk[uint8(i-1)]
            sk[pixel] = temp + prob[pixel]
        }
    }
    return sk
}

The entire pixel range from 0-255 will be stored in the map as keys, the value of each key is its cdf. So there will be values from 0 to ~1. The cdf function above is continued by multiplying each value by its maximum value. This is done to convert the CDF into a full range of values, i.e. from 0 to the desired maximum value. In this way, it is expected to expand the pixel intensity distribution so that it covers the entire possible range.

func (ip *ImageProcessor) HistogramEqualizationLevel(cdfMultiply map[uint8]float64) map[uint8]int {
    histEQ := make(map[uint8]int)

    for i := 0; i <= 255; i++ {
        histEQ[uint8(i)] = int(cdfMultiply[uint8(i)] + 0.5)
    }

    return histEQ
}

The function above handle finding the desired histogram equalization level. The number 255 is a static number that I use because the maximum number of 8bit integers is 255, or values from 0-255. That's it.

The next step is the last step.

func (ip *ImageProcessor) ApplyHistogramEqualization(originalImage [][]color.Gray, eqLevel map[uint8]int) [][]color.Gray {
    height := len(originalImage)
    width := len(originalImage[0])

    equalizedImage := make([][]color.Gray, height)
    for y := 0; y < height; y++ {
        equalizedImage[y] = make([]color.Gray, width)
        copy(equalizedImage[y], originalImage[y])
    }

    for y := 0; y < height; y++ {
        for x := 0; x < width; x++ {
            pixel := originalImage[y][x]
            equalizedPixelValue := eqLevel[uint8(pixel.Y)]
            equalizedImage[y][x].Y = uint8(equalizedPixelValue)
        }
    }

    return equalizedImage
}

Finally, there is the process of applying histogram equalization to the grayscale image. The first thing to do is to take the dimensions, followed by building an image with predetermined dimensions. The last process is the application of histogram equalization to the previously created image.

At this point, the histogram equalization process should have been completed. But in order to be formed as a.Gray image on the golang, I created a helper function for the equalized data to be applied to the grayscale image.

func (ip *ImageProcessor) CreateImageFromGrayPixels(pixels [][]color.Gray) *image.Gray {
    height := len(pixels)
    width := len(pixels[0])

    img := image.NewGray(image.Rect(0, 0, width, height))

    for y := 0; y < height; y++ {
        for x := 0; x < width; x++ {
            img.Set(x, y, pixels[y][x])
        }
    }

    return img
}

Final output of the histogram equalization implementation below.

It can be seen in the output that the image appears sharper based on what the eye sees. The mathematical proof is that, in the histogram, the range changes from the original approximately 50-200 pixel variation, after equalization the histogram becomes 0 to 250, indicating that the pixel intensity distribution has been changed such that the contrast between the various pixel intensity values is enhanced. By increasing the contrast, the difference between dark and light areas in the image becomes more pronounced, so the image looks sharper and the details are better defined. This results in a more distinct difference between objects in the image, which makes it more recognizable to the human eye.

You can see the full of code on github

Intepretasi Regresi Logistik: Dari Logit Hingga Penurunan Gradient.

Rama Reksotinoyo — Sun, 08 Oct 2023 07:25:46 +0000

Dengan menyebut yang dipertuan agung programmer andal di luar sana dan ChatGPT, saya haturkan salam sebelum menuliskan kode program pada artikel ini.

Regresi logistik merupakan satu dari banyak algoritma machine learning yang sering dipakai untuk kasus klasifikasi. Intepretasi model regresi logistik pada curva membentuk melengkung/non-linear. Ini adalah kurva regresi logistik.

Model diatas dapat diintepretasikan sebagai garis yang memilki kemiringan untuk βπ(x)[1-π(x)] terhadap parameter β.

Jika persamaan diatas dielaborasi lebih dalam lagi, π(x) menggambarkan probabilitas suatu peristiwa terjadi pada titik x. Persamaan ini mengasumsikan bahwa probabilitas peristiwa tersebut dapat diprediksi atau diperkirakan menggunakan fungsi π(x). Misalnya, jika π(x) memiliki nilai 0,5, maka probabilitas peristiwa terjadi pada titik x adalah 0,5. Dan titik 0,5 sendiri merupakan titik yang paling curam pada kurva.

Masuk ke bagian yang paling penting, yaitu intepretasi dari model yang dihasilkan. Hasilnya akan direpresentasikan dalam bentuk logit atau log odds, dengan odds sendiri adalah peluang sukses dibagi dengan peluang gagal. Dengan persamaan dibawah ini:

log(odds) = β₀ + β₁x
Intepretasi : setiap unit pertambahan pada variabel x akan mengakibatkan pertambahan pada β1 unit pada log odds.

Selanjutnya, untuk mengintepretasikan oddds, maka akan ditransformasi nilai log dari odds menjadi odds.
Dengan demikian,

odds = e^(β₀ + β₁x) = e^(β₀) * (e^(β₁))^x

Interpretasi: odds perkiraan keberhasilan akan dikalikan dengan exp(β₁) untuk setiap peningkatan satu unit pada x.

Masuk kedalam kode program, disini saya mencoba mengimplementasikan regresi logistik dengan menggunakan bahasa pemrograman C yang tentu saja dibantu oleh ChatGPT. Dataset yang saya gunakan adalah dataset mengenai horseshoe crab atau masyarat Indonesia lebih mengenal dengan sebutan belangkas. Dataset berisi satu fitur atau prediktor berupa lebar cangkang, dan target kelas atau respon berisi punya atau tidaknya satelite atau tanduk pada kepiting ini. Jumlah record pada dataset adalah 173 baris dengan data sudah bersih dari data kosong, dublikat, dan ketidak konsistenan data.

Pertama yang saya lakukan adalah membuat struct yang nantinya akan direpresentasikan sebagai objek, baik itu sebagai objek data atau pun ndarray (N-Dimensional Array) atau array berdimensi tinggi.

typedef struct {
  unsigned int num_rows;
  unsigned int num_cols;
  float **data;
  int is_square;
} Ndarray;

Lalu saya membuat fungsi yang digunakan untuk memuat dataset yang nantinya digunakan pada program ini.

Ndarray* readCSV(const char* filename, int maxRecords)
{
  FILE* file = fopen(filename, "r");
  if (file == NULL)
  {
    printf("Error opening file.\n");
    return NULL;
  }

  int records = 0;
  int read = 0;
  float** data = malloc(maxRecords * sizeof(float*));
  if (data == NULL) {
    printf("Error allocating memory.\n");
    fclose(file);
    return NULL;
  }

  while (records < maxRecords)
  {
    data[records] = malloc(2 * sizeof(float));
    if (data[records] == NULL) {
      printf("Error allocating memory.\n");
      for (int i = 0; i < records; i++) {
        free(data[i]);
      }
      free(data);
      fclose(file);
      return NULL;
    }

    read = fscanf(file, "%f,%f\n", &data[records][0], &data[records][1]);

    if (read == 2)
      records++;
    else if (read != 2 && !feof(file))
    {
      printf("File format incorrect.\n");
      for (int i = 0; i < records; i++) {
        free(data[i]);
      }
      free(data);
      fclose(file);
      return NULL;
    }

    if (feof(file))
      break;
  }

  fclose(file);

  Ndarray* crabs = malloc(sizeof(Ndarray));
  if (crabs == NULL) {
    printf("Error allocating memory.\n");
    for (int i = 0; i < records; i++) {
      free(data[i]);
    }
    free(data);
    return NULL;
  }

  crabs->num_rows = records;
  crabs->num_cols = 2;
  crabs->data = data;

  return crabs;
}

Demi mengurangi drama type casting, saya terpaksa membuat varibel yang berupa bilangan real saya rubah ke desimal.

Dilanjutkan dengan saya membuat fungsi untuk mengetahui berapa jumlah sukses(y=1) dan jumlah gagal pada dataset(y=0).

int get_n_success(Ndarray* crabs)
{
  int n_success = 0;

  for (int i = 0; i < crabs->num_rows; i++)
  {
    if (crabs->data[i][0] == 1.0)  // Perubahan pada baris ini
    {
      n_success += 1;
    }
  }

  return n_success;
}

int get_n_failed(Ndarray* crabs)
{
  int n_failed = 0;

  for (int i = 0; i < crabs->num_rows; i++)
  {
    if (crabs->data[i][0] == 0.0)
    {
      n_failed += 1;
    }
  }

  return n_failed;
}

Hasil dari kode program diatas adalah n_success sebanyak 111, dan n_failure sebanyak 62.

Selanjutnya adalah gradient decent, yang saya yakini dengan subjektifitassaya merupakan salah satu algoritma terbaik yang pernah diciptakan manusia. Yang mana, fungsi ini akan melakukan optimisasi parameter, pada kasus ini parameter hanya lebar cangkang.

Karena tujuan dari gradient descent adalah memperbarui paramenter, maka untuk mencari nilai optimum dari β₀ and β₁, memaksimalkan fungsi likelihood/log-likelihood, atau meminimalkan negative log-likelihood:
$log(L(β)) = -Σ{i=1}^{n}(y{i}log[π(x_{i})]+(1-y_{i})log[1-π(x_{i})])$
Dan jika parameternya sebanya j, maka:

Untuk mencari parameter βj atau parameter sebanyak j maka:

Selajutnya sebelum membangun fungsi gradient decent, perlu disiapkan juga fungsi non-linear sigmoid, untuk mengintepretasi hasil klasifikasi, entah hasilnya akan masuk ke kelas 0 atau 1. Dibawah ini adalah implementasi dari fungsi sigmoid:

double sigmoid(double *X, int num_cols, double b0, double *b1) {
    double logit = b0;
    double exp_val, pi;

    for (int i = 0; i < num_cols; i++) {
        logit += X[i] * b1[i];
    }

    exp_val = exp(logit);
    pi = exp_val / (1 + exp_val);

    return pi;
}

Dilanjutkan dengan membangun fungsi yang dapat menghitung nilai sigmoid disetiap data point:

Ndarray* calculate_sigmoid_list(Ndarray* matrix, double b0, double *b1) {
    Ndarray* sigmoid_list = malloc(sizeof(Ndarray));
    sigmoid_list->num_rows = matrix->num_rows;
    sigmoid_list->num_cols = 1;
    sigmoid_list->data = malloc(sigmoid_list->num_rows * sizeof(double*));

    for (int i = 0; i < sigmoid_list->num_rows; i++) {
        sigmoid_list->data[i] = malloc(sigmoid_list->num_cols * sizeof(double));

        double* matrix_row = (double*) matrix->data[i];

        sigmoid_list->data[i][0] = sigmoid(matrix_row, matrix->num_cols, b0, b1);
    }

    return sigmoid_list;
}

Langkah selanjutnya adalah membuat fungsi kerugian atau nilai error, pada impelemtasinya, saya menerapkan negative log loss-likelihood. Dengan persamaan:
$log(L(β)) = -Σ{i=1}^{n}(y{i}log[π(x_{i})]+(1-y_{i})log[1-π(x_{i})])$
Dan fungsi pada kode programnya adalah:

double cost_function(Ndarray* matrix, double *pi) {
    double eps = 1e-10;
    double log_loss = 0.0;

    for (int i = 0; i < matrix->num_rows; i++) {
        double log_like_success = matrix->data[i][0] * log(pi[i] + eps);
        double log_like_failure = (1 - matrix->data[i][0]) * log(1 - pi[i] + eps);
        double log_like_total = log_like_success + log_like_failure;
        log_loss -= log_like_total;
    }

    return log_loss;
}

Hasilnya adalah saya mendapat nilai kerugian sebesar:
119.91446220227053

Lalu saya menyiapkan fungsi untuk menghitung turunan dari β₀ atau intercept. Dengan persamaan:
$gradβ₀=(∑i=1n[π(xi)−yi])𝑥1grad_β₀ = (\sum_{i=1}^{n} [\pi(x_{i}) - y_{i}])𝑥1$

Kode programnya adalah:

float gradient_b0(Ndarray* matrix, double *pi){
  float grad;
  for(int i=0; i<matrix->num_rows; i++){
    grad += (pi[i] - matrix->data[i][0]);
  }

  return grad;
}

Sedangkan persamaan matematis dan fungsi untuk menghitung gradient dari β₁ adalah:
$gradβ₁=(∑i=1n[π(xi)−yi])𝑥grad_β₁ = (\sum_{i=1}^{n} [\pi(x_{i}) - y_{i}])𝑥$

double* gradient_b1(Ndarray* matrix, double* pi, double* y, double b0, double* b1) {
    int num_rows = matrix->num_rows;
    int num_cols = matrix->num_cols;
    double* grad_b1 = (double*)malloc(num_cols * sizeof(double));

    for (int j = 0; j < num_cols; j++) {
        grad_b1[j] = 0.0;
        for (int i = 0; i < num_rows; i++) {
            double* matrix_row = (double*)matrix->data[i];
            double sigmoid_val = sigmoid(matrix_row, num_cols, b0, b1);
            grad_b1[j] += matrix_row[j] * (sigmoid_val - y[i]);
        }
    }

    return grad_b1;
}

Setelah fungsi menghitung gradient dari β₀ dan β₁ dibuat, langkah terakhir dalam gradient descent adalah menerapkan langkah-langkah dibawah ini:

Algoritma Optimasi:

Definisikan variable response y and variabel prediktor X.
Inisialisasi parameter yang akan optimasi β₀=0.0 and β₁=0.0.
Hitung gradient dari ∂NLL/∂𝛽𝑗, sehingga mendapatkan gradient dari β₀ and β₁.
Perbarui paremeter β₀ and β₁ $βj^{new} = βj^{old} - lr *$ ∂NLL/∂𝛽𝑜𝑙𝑑𝑗
Ulangi langkah 2-4: Δ𝛽𝑗 < nilai yang ditoleransi, or ∇𝛽𝑗NLL < nilai yang ditoleransi.

Kode program dari langkah-langka diatas adalah:

void gradient_descent(double *X, double *y, int num_rows, int num_cols, double eta, double tol, double *b0, double *b1, int *iterations, double *log_loss_list) {
    srand(time(NULL));
    double b1_initial[2] = {((double)rand() / RAND_MAX) * 0.01, ((double)rand() / RAND_MAX) * 0.01};

    // Make a criteria to run the iteration
    int continue_iteration = 1;

    // Running the iteration
    int i = 0;
    while (continue_iteration) {
        // Update i
        i++;

        // Calculate success probability (pi) from the current b0 and b1
        Ndarray matrix;
        matrix.num_rows = num_rows;
        matrix.num_cols = num_cols;
        matrix.data = malloc(matrix.num_rows * sizeof(double*));
        for (int j = 0; j < matrix.num_rows; j++) {
            matrix.data[j] = malloc(matrix.num_cols * sizeof(double));
            for (int k = 0; k < matrix.num_cols; k++) {
                matrix.data[j][k] = X[j * matrix.num_cols + k];
            }
        }
        Ndarray* sigmoid_list = calculate_sigmoid_list(&matrix, *b0, b1_initial);
        double *pi = malloc(num_rows * sizeof(double));
        for (int j = 0; j < num_rows; j++) {
            pi[j] = sigmoid_list->data[j][0];
        }

        // Calculate log loss from the current b0 and b1
        double log_loss = cost_function(&matrix, pi);
        log_loss_list[i - 1] = log_loss;

        // Calculate gradient of b0 and b1
        double grad_b0 = gradient_b0(&matrix, pi);
        double* grad_b1 = gradient_b1(&matrix, pi, y, *b0, b1_initial);

        // Update b0 and b1
        *b0 -= eta * grad_b0;
        for (int j = 0; j < num_cols; j++) {
            b1[j] -= eta * grad_b1[j];
        }

        // Check convergence
        if (i > 1 && fabs(log_loss_list[i - 1] - log_loss_list[i - 2]) < tol) {
            continue_iteration = 0;
        }

        // Free memory
        free(pi);
        free(grad_b1);
        free(sigmoid_list->data);
        free(sigmoid_list);
        for (int j = 0; j < matrix.num_rows; j++) {
            free(matrix.data[j]);
        }
        free(matrix.data);
    }

    // Store the number of iterations
    *iterations = i;

    // Print final b0 and b1
    printf("Nilai akhir b0: %f\n", *b0);
    for (int j = 0; j < num_cols; j++) {
        printf("Nilai akhir b1[%d]: %f\n", j, b1[j]);
    }
}

Hasil dari nilai optimum dari β₀ dan β₁ berturut-turut adalah:
-7.971823 dan 0.009542.

Sehingga logit model dengan parameter yang optimum adalah:
log(odds(𝜋(𝑥))=-7.971823+0.009542(𝑥)

Kesimpulan dari model logit diatas adalah:

Odds Ratio: Koefisien regresi logistik (0.009542 dalam persamaan ini) mengukur perubahan dalam log-odds dari variabel prediktor (𝑥) terhadap kemungkinan kejadian sukses (atau kategori positif) dalam masalah klasifikasi biner. Dalam hal ini, setiap peningkatan satu unit dalam 𝑥 akan menghasilkan peningkatan sebesar 0.009542 dalam log-odds dari kejadian sukses. Dalam kasus Anda, variabel 𝑥 mungkin adalah salah satu fitur atau prediktor yang digunakan dalam model.
Intercept (Konstanta): Intercept (-7.971823 dalam persamaan ini) adalah nilai log-odds ketika semua variabel prediktor lainnya adalah nol. Dalam interpretasi praktis, ini mewakili dasar atau tingkat dasar log-odds kejadian sukses ketika tidak ada pengaruh dari variabel prediktor lainnya.
Log-Odds: Persamaan ini menggambarkan logaritma dari odds (rasio peluang) bahwa kejadian sukses terjadi. Dalam beberapa kasus, Anda mungkin perlu mengonversi log-odds menjadi peluang atau probabilitas jika Anda ingin mendapatkan interpretasi yang lebih intuitif.
Pengaruh Prediktor: Koefisien regresi (0.009542) mengindikasikan seberapa kuat pengaruh 𝑥 terhadap peluang kejadian sukses. Jika nilai koefisien positif, maka peningkatan dalam 𝑥 akan meningkatkan log-odds kejadian sukses, sementara jika nilai koefisien negatif, maka peningkatan dalam 𝑥 akan mengurangi log-odds kejadian sukses.

Dapat dihitung nilai odds ratio e^0.009542 untuk mengetahui seberapa besar efek 𝑥 terhadap peluang kejadian sukses. Jika odds ratio lebih besar dari 1, maka peningkatan dalam 𝑥 akan meningkatkan peluang kejadian sukses. Jika odds ratio kurang dari 1, maka peningkatan dalam 𝑥 akan mengurangi peluang kejadian sukses.

Kode program lengkap saya tulis pada github

Saran pengembangan selanjutnya:

Visualisasi hasil gradient descent.
Menghitung interval probability.
Melakukan significance testing terhadap kemungkinan seekor kepiting memiliki satelit tidak bergantung pada lebarnya atau predictornya.
Melakukan standarisasi pada fitur atau predictor.

Referensi:
Categorical Statistic Course Pacmann

Writing your own linear algebra matrix library in C

An Introduction to Statistical Learning

AB Testing Ads - Project Experimental & AB Testing

Rama Reksotinoyo — Wed, 06 Sep 2023 13:56:41 +0000

Perusahaan pemasaran ingin menjalankan kampanye yang sukses, namun pasarnya kompleks dan dengan pilihan kampanye yang tepat, mereka dapat menghasilkan kampanye yang sukses. Jadi mereka biasanya melakukan pengujian A/B. Kebanyakan orang akan dihadapkan pada iklan (kelompok uji). Dan sejumlah kecil orang (kelompok kontrol) akan melihat iklan layanan masyarakat (PSA) dengan ukuran dan lokasi yang sama dengan iklan tersebut.

Data yang digunakan pada eksperimen ini didapat dari kaggle marketing AB Testing. Fitur dari dataset terdiri dari :

User ID (unique)
Test group : Jika “ad” orangnya melihat iklan, kalau “psa” hanya melihat iklan.
Converted : Jika seseorang membeli produk true, jika tidak else.
Total ads : Jumlah iklan.
Most ads day : Hari dimana seseorang melihat iklan terbanyak.
Most ads hour : Jam dimana seseorang melihat iklan terbanyak.

Tujuan eksperimen yang dilakukan adalah untuk menentukan apakah iklan pemasaran sukses atau tidak. Dengan mengetahui hal tersebut, diharapkan kedepannya perusahaan dapat mengoptimalkan pemasaran.

Metriks yang digunakan pada eksperimen ini adalah conversion rate. Dengan group control public service announcement dan group treatment akan ada iklan dengan warna yang mencolok.

Hipotesis yang akan diuji pada eksperimen ini adalah :

H0: tidak terdapat perbedaan yang signifikan dari implementasi ads terhadap conversion rate.
H1: terdapat perbedaan yang signifikan antara grup control dan treatment akibat dari implementasi ads.

Untuk bagian desain eksperiment, randomization unitnya adalah UserID karena unik, target of randomization unit adalah user atau pengguna secara acak. Pada sample size significant level atau alfa yang telah ditetapkan pada eksperiment kali ini adalah 5% atau 0.05. Power level atau probabilitas untuk menolak hipotesis null, yang menggunakan nilai 0.08. Standar deviasi dari populasinya adalah 0.07. Yang terakhir adalah Difference between control and treatment atau delta adalah 1% atau 0.01.

Untuk kode program, pertama yang dilakukan adalah memeriksa adakah data yang kosong atau duplikat, setelah dicek tidak ada data yang kosong ataupun duplikat.

Selanjutnya, dilakukan perhitungan untuk mencari jumlah ukuran sample, dengan cara sebagai berikut:

def sample_size_calculator(alpha: float, power_level: float, standard_deviation: float, delta: float) -> int:
    beta = 1 - power_level
    z_alpha = stats.norm.ppf(1 - alpha/2)
    z_beta = stats.norm.ppf(1 - beta)
    n = (2 * (z_alpha + z_beta)**2 * standard_deviation**2) / delta**2
    return int(n)

Dari kode program diatas didapat hasil ada 769.

Untuk durasi eksperimen, dengan mempertimbangkan jumlah pengunjung yang sekiranya dalam 2 minggu terakhir akan terus menggunakan website, maka eksperimen akan berlangsung sekitar 2 minggu.

Tahap selanjutnya yang dilakukan adalah dengan membagi dataset utk grup treatment dan control, hasilnya masing-masing grup treatment dan control sebanya 769.

data_control = df[df['test_group'] == 'psa'].sample(sampe_size)
data_treatment = df[df['test_group'] == 'ad'].sample(sampe_size)

len(data_control), len(data_treatment) # (769, 769)

Lalu saya melakukan visualisasi untuk melihat distribusi total ads dan most ads user

Selanjutnya adalah mencari tahu jumlah user yang converted baik pada control maupun treatment:

n_user_control = data[data["test_group"] == "psa"]["user_id"].nunique()
n_user_treatment = data[data["test_group"] == "ad"]["user_id"].nunique()
# 769

n_converted_control = data[(data["test_group"] == "psa") & (data["converted"] == True)].shape[0]
n_converted_treatment = data[(data["test_group"] == "ad") & (data["converted"] == True)].shape[0]

Lalu menghitung nilai observed dan expected-nya:

# Menghitung total jumlah pengguna yang terlibat
total_user_impacted = n_user_control + n_user_treatment

# Menghitung nilai observed dan expected
observed = [n_user_control, n_user_treatment]
expected = [int(total_user_impacted / 2), int(total_user_impacted / 2)]

Lalu mencari convert rate untuk control dan treatment:

cr_control = n_converted_control / n_user_control
cr_treatment = n_converted_treatment / n_user_treatment

Langkah selanjutnya adalah uji hipotesis, utk itu saya menyiapkan variable variable yang diperlukan seperti z_score, z_crit, dan juga p_value.

Hasilnya adalah sebagai berikut

Pada gambar diatas H0 ditolak, karena group treatment lebih besar dari pada group control, yang mana hal ini berakibat bahwa treatment berpengaruh pada pemasaran iklan.

Pada uji hipotesisi, terlihat bahwa adanya perbedaan yang terlihat antara grup treatment dan grup control, sehingga iklan diterima dengan baik oleh sebagian pelanggan pada kelompok treatment.

Reference:
Garson, G.D., 2012. Testing statistical assumptions: Blue Book Series. Asheboro: Statistical Associate Publishing, [online] pp.1–54

The Divine Distribution: A Hands-on, Humorous and Spiritual Take on the Central Limit Theorem.

Rama Reksotinoyo — Sat, 22 Jul 2023 19:32:34 +0000

There is an ascetic from Tibet who says that "If several prerequisites are satisfied, the distribution of the mean of multiple independent random variables tends to approximate a normal distribution as the size of the sample tends to be large."

I joked, that was the sound of the central limit theorm that was hard to digest, I got it from the GPTModel :).

The formal definition of the central limit theorem, as is commonly heard in classrooms, is

"The Central Limit Theorem is a statistical principle that claims that no matter how the population was originally distributed, if we take numerous random samples from it, the distribution of the sample means will resemble a normal distribution. In other words, the closer the distribution of the sample mean is to the normal distribution, the more samples are taken."

One of the applications of the quotation above when doing sampling, where when I do sample with the terms of each sample data > 30 then the normal tendency distribution will be seen. No matter what the distribution on its population, whether it is skewed positive or negative, the result will remain the same, that is, will be in the form of normal distribution.

I will start creating a data distribution with the amount of data as much as 10000 times, with its descriptive statistics being an average of 90, the spread of data seen with the devisation standard of 10. I implemented it bellow:

Before going too far, I implemented the content in this article using the Python programming language, with my library I prepared:

numpy
pandas
matplotlib
seaborn
scipy

The visualization of the distribution:

After that, lets try the theorm, first i created a function for generate a sample. I tried to sample 10,000 times, with each sample containing 25 data, resulting in an average of 89,981172 and a standard deviation of 1,999175. Here's the function:

The mathematical equation for calculating the sample mean is as follows:

Even when plotting is done, the sample that I have generated will produce a distribution that is more or less the same as the population or normally distributed.

Using the generate_population_data function, I add the 180 degree parameters of skewnes and also positive type skewness, and then I move on to the fascinating part, which is trying to do sampling on populations that are distributed skewed whether it is positive or negative. The distribution visualization is as follows:

From the sample distribution above, I sampled using the generate_sample function, by sampling 3000 times and each sample has 30 data. The average result is 11.103568 and the standard deviation is 10.807124. Maybe the output generated by the function will be different every time you run the program code, but it's unlikely to be very different.
The visualization of the sampling distribution is as follows:

here's the generate_sample function.

Following successively is a visualization of each population distribution generated by the generate_population_data function and a visualization of each sample taken from each population distribution using generate_sample function.

The Divine Distribution, a concept with a humorous and spiritual twist, is proven to be a truth in statistics. When numerous independent random variables are sampled and the sample size is sufficiently large while fulfilling certain conditions, the distribution of the mean tends to follow a normal distribution. This aligns with the divine nature of the normal distribution, which is perfect and can be applied universally, regardless of the shape of the original population distribution. Therefore, the Divine Distribution is a powerful tool that enables us to make accurate predictions and inferences about a population based on samples.

The code is available on github.

References:

Pacmann course
Wikipedia
Statquest

Analisis Data Tagihan Kesehatan - Project Probability Pacmann.

Rama Reksotinoyo — Sun, 09 Jul 2023 07:21:20 +0000

Pada project akhir kelas probability Pacmann, saya melakukan analisis terhahadap beberapa variabel. Project ini bertujuan untuk mengetahui hubungan kesehatan dengan biaya kesehatan.

Pertama dataset yang digunakan memiliki field yang berisi seperti dibawah ini:

age: Usia
bmi: Indeks masa tubuh
children: Anak yang ditanggung oleh asuransi kesehatan
smoker: Merokok atau tidaknya
region: Area tempat tinggal
charges: Biaya yang yang ditagih oleh asuransi kesehatan

Beberapa outline yang saya kerjakan diantaranya adalah menampilkan statistik deskriptif pada dataset, pmf, pdf, korelasi antar dua variabel, dan hypotesis testing.

Pada statistik deskriptif saya menampilkan beberapa rata-rata dari beberapa kolom, seperti dibawah ini:

Selanjunya adalah analisa peluang pada variable diskrit. Beberapa peluang yang saya telah analisis adalah:

Selanjutnya saya lanjutkan dengan menganalisa variabel continue, yang mana untuk mengethui peluang "Seseorang dengan BMI di atas 25 mendapatkan tagihan kesehatan di atas 16.7k, atau Seseorang dengan BMI di bawah 25 mendapatkan tagihan kesehatan di atas 16.7k"

Selanjutnya, saya menganalisa korelasi variabel sex, smoker, dan region. Langkahnya adalah dengan menampilakn visualisasi korelasi menggunakan heatmap dan dilanjutkan dengan membandingkan nilai korelasi antara variabel age dengan charges, bmi dengan charges, dan children dengan charges.

Terakhir yang saya bahas pada project akhir probability saya adalah mengenai hypotesis testing, atau cara untuk menguji kevalidan suatu hypotesis. Beberapa yang saya uji adalah:

Tagihan kesehatan perokok lebih tinggi daripada tagihan kesehatan non perokok
Tagihan kesehatan laki-laki lebih besar dari perempuan
Proporsi perokok laki laki lebih besar dari perempuan

github

references:

Materi kelas pacmann

Analisis data pada olist database - Project Data Wrangling Pacmann.

Rama Reksotinoyo — Thu, 06 Jul 2023 16:21:55 +0000

Hasil dari yang telah saya baca secara singkat di internet, database yang digunakan pada projct ini adalah database pada suatu toko online, yang mencakup data barang, transaksi, hingga penjual. Olist database terdiri dari 8 buah tabel, seperti yang telah saya lihat menggunakan tools pandas.

Tetapi saya tidak akan menggunakan keseluruh tabel diatas, hanya menggunakan beberapa tabel yang akan digunakan sesuai dengan objective yang telah saya tentukan pada awal pengerjaan proyek.

Objective yang saya tentukan adalah:

Jika suatu perusahaan ingin mengetahui 10 produk yang paling laku.
Jika perusahaan ingin mengetahui 10 produk yang banyak di-cancel.

Tahap pertama adalah saya melihat isi dari tabel yang saya butuhkan dengan menggunakan "SELECT * FROM ....". Setelah itu saya melakukan preprocessing. Tahap preprocessing yang saya lakukan adalah mengganti product name bahasa spanyol menjadi bahasa inggris dengan menggunakan left join dengan tabel category products. Dilanjutkan dengan melakukan mengganti manual value produk yang tidak ada di tabel product_category. Prosesnya seperti dibawah ini:

Proses pemeriksaan selanjutnya adalah dengan mengecek ada atau tidaknya missing value pada data. Kolom yang ada missing value-nya memilki perlukan berbeda, tergantung konteks. Pada konteks project kali ini mengisi kolom kosong dengan pada tabel produk, seperti pada tabel dibawah ini.

Selanjutnya adalah proses pengecekkan ada atau tidaknya baris yang duplikat antara data satu dengan yang lain, sama atau tidaknya diukur dari baris yang berurutan apakah sama persis value antara kolomnya atau tidak. Dilakukan pengecekan duplicate karena banyak ditemui data-data yang seharusnya duplicate, tetapi dibiarkan saja.

Tahap selanjutnya adalah melakukan pengecekkan data dari kolom product_name dari tabel new_product. Prosesnya adalah seperti gambar dibawah ini

Proses selanjutnya adalah dilakukannya pengecekan ada atau tidaknya outlier pada kolom harga. Visualisasi dilakukan dengan boxplot agar data outlier dan data-data yang lain bisa dengan mudah terlihat.

Karena banyaknya outlier, maka saya rubah ke dalam bentuk log. Hasil dari visualisasi menggunakan boxpot adalah seperti ini.

Selanjutnya adalah proses join table agar nantinya pada saat analisis tidak perlu melakukan join yang akan menyulitkan pandangan pembaca.

Dan yang terakhir adalah melihat objektif yang telah saya tentukan diawal, pada proses ini dilakukan visualisasi dengan barplot dengan menggunakan library seaborn . Diantarnya adalah sebagai berikut:

Produk yang paling banyak dijual (Menampilkan 10 produk).

Mengetahui produk apa saya yang banyak di-cancel oleh pelanggan.

github

Walking through the bayes theorm to update our belief systematically.

Rama Reksotinoyo — Sun, 28 May 2023 17:29:20 +0000

I have heard about a quote from ancient Greek that goes, and this quote bellow is related to Bayes' theorem.

"Always leave room for uncertainty so that we remain open to new things and survive."

Bayes' theorem is a statistical method used to calculate the uncertainty or probability of an event based on past assumptions that may be related to the event. Inspired by this, I decided to delve deeper (with my limited brain capacity, of course) into Bayes' theorem, which I tried to implement in program code.

I attempted to implement the probability of someone, let's call them X, having HIV and actually being HIV positive. I wanted to determine the probability that X truly has HIV.

Here is the formula for Bayes' theorem:

With my zero of medical knowledge, I assumed that out of a population of 5000 people who took an HIV test, only 3% of the population actually had HIV, considering that I rarely come across people with HIV in my current environment.

After reading some papers and hearing opinions from HIV specialists and HIV researchers, I found that if X were to take the test now, the probability of them having HIV would be 95%. This can be referred to as sensitivity. If X does not have HIV, the probability of them truly not having HIV is 98%. This can be referred to as specificity.

Therefore:
P(HIV) = prior = 0.3
P(diagnosed|HIV) = likelihood = 0.95
P(diagnosed) = evidence
P(HIV|diagnosed) = posterior
P(not HIV|not diagnosed) = 0.98

Thus, the confusion matrix table is as follows, where the horizontal axis represents the actual values and the vertical axis represents the predicted values.

Please note that the translation provided may not precisely match the terminology used in statistical literature, but it conveys the general meaning.

Confusion Matrix:

Since the evidence is unknown in this case, I am trying to calculate it using the equation from the law of total probability, which is P(B) = Σ[P(B|Aᵢ)P(Aᵢ)]. Here, I am implementing it using the Go programming language:

const (
    PDiagnosedGivenHIV      = 0.95 // P(D|S)
    PNotDiagnosedGivenNoHIV = 0.98 // P(-D|-S)
    PDiagnosedGivenNoHIV    = 0.02 // P(D|¬S)
    PNotDetectedGivenSpam   = 0.05 // P(-D|S)
)

func TotalProbability(prior *float64) float64 {
    // using law of total probability.
    // P(D) = P(D|S)P(S) + P(D|¬S)P(¬S)
    PNoHIV := 1 - *prior
    return (PDiagnosedGivenHIV * *prior) + (PDiagnosedGivenNoHIV * PNoHIV)
}

With the obtained evidence of 0.0479, I just need to substitute the above evidence with the known data into the following equation.

Here is the complete code along with the test cases:

package main

import "fmt"

type Probability struct {
    PriorProb       float64
    ConditionalProb float64
    EvidenceProb    float64
}

func (p *Probability) CalculatePosteriorProb() float64 {
    // Menghitung posterior probability menggunakan rumus Bayes
    posteriorProb := (p.ConditionalProb * p.PriorProb) / p.EvidenceProb
    return posteriorProb
}

const (
    PDiagnosedGivenHIV      = 0.95 // P(D|S)
    PNotDiagnosedGivenNoHIV = 0.98 // P(-D|-S)
    PDiagnosedGivenNoHIV    = 0.02 // P(D|¬S)
    PNotDetectedGivenSpam   = 0.05 // P(-D|S)
)

func TotalProbability(prior *float64) float64 {
    // using law of total probability.
    // P(D) = P(D|S)P(S) + P(D|¬S)P(¬S)
    PNoHIV := 1 - *prior
    return (PDiagnosedGivenHIV * *prior) + (PDiagnosedGivenNoHIV * PNoHIV)
}

func main() {
    prior := 0.03
    evidence := TotalProbability(&prior)

    prob := Probability{
        PriorProb:       prior,    // Probabilitas prior
        ConditionalProb: 0.95,     // Probabilitas kondisional
        EvidenceProb:    evidence, // Probabilitas evidence
    }

    // Menghitung posterior probability menggunakan Bayes' Rule
    posteriorProb := prob.CalculatePosteriorProb()

    fmt.Printf("Posterior Probability: %.2f\n", posteriorProb)
    fmt.Println(" ")

}

package main

import (
    "math"
    "testing"
)

const epsilon = 0.00000001 // Toleransi kesalahan

func TestCalculatePosteriorProb(t *testing.T) {
    prob := Probability{
        PriorProb:       0.1,  // Probabilitas prior
        ConditionalProb: 0.75, // Probabilitas kondisional
        EvidenceProb:    0.1,  // Probabilitas evidence
    }

    expectedPosteriorProb := 0.75

    actualPosteriorProb := prob.CalculatePosteriorProb()

    if math.Abs(actualPosteriorProb-expectedPosteriorProb) > epsilon {
        t.Errorf("Wrong answer! it should be %.8f", expectedPosteriorProb)
    }
}

func TestTotalProbability(t *testing.T) {

    const (
        PDetectedGivenSpam       = 0.95 // P(D|S)
        PNotDetectedGivenNonSpam = 0.98 // P(-D|-S)
        PDetectedGivenNonSpam    = 0.02 // P(D|¬S)
        PNotDetectedGivenSpam    = 0.05 // P(-D|S)
    )

    e := 0.0001
    prior := 0.2

    expectedTotalProb := 0.2060
    actualTotalProb := TotalProbability(&prior)

    if math.Abs(actualTotalProb-expectedTotalProb) > e {
        t.Errorf("Wrong answer! It should be %.4f but %.4f", expectedTotalProb, actualTotalProb)
    }
}

In the next articel , I will talk about how to update our level of confidence or prior, which will demonstrate that the likelihood of someone being affected by HIV is directly proportional to the magnitude of our prior or initial assumption. Testing and implementation of the narrative "to update our belief systematically" will be in the next section. See you soon.

Reference:
Wikipedia

Zenius

StatQuest

Automatic Differentiation with Finite Difference Approximation. That’s what it is.

Rama Reksotinoyo — Thu, 15 Dec 2022 18:14:48 +0000

Here is a gentle introduction about automatic differentiation in modern machine learning libraries or deep learning libraries like Pytorch, Jax, or even tensorflow. One of the powerful features of them is auto gradient, or we don't have to initiate the function blablabla with respect to blablabla, because the machine automatically computes the gradient automatically. The implementation of automatic gradient can be seen when we compute the gradient descent to minimize the error with respect to parameters.

So, what exactly the auto differential is? if you jump into this articles without machine learning background you may be confusing about auto diff or auto gradient, in this article I will consider it the same, and for now I will use auto gradient terms. I’ll repeat it, what exactly auto gradient is? basically auto gradient is a numerical method for efficiently calculating the derivative of any function using programming techniques.

For this part I will reduce the mathematical terms or explanations or whatever, cause i'm not a mathematical guy, neither knows the mathematical idea behind deeply. So I will explain as a self-proclaimed ordinary engineer

Came up from the definition of the derivative above, which can then be elaborated again with the equation.

If the value of the xi element is slightly added and subtracted while the value of the other elements stays the same, the quantifier on the right side represents the difference in the value of the function f. The finite difference method is the name of the technique used to calculate the derivative's value.

And also, we talked about derivatives, whether native derivative, partial derivative, and chain rule derivative. Those three kinds of derivatives are needed for this idea.

It should be mentioned that Larry Gonick discusses derivatives in a very amiable manner in his book Cartoon Calculus.

In machine learning or deep learning, auto gradients can be used for optimization. Or we use the best algorithm ever, gradient descent (obviously, it's a hot take for me). With the optimization, machine learning model can learn and the updates the parameters while the process of back-forward propagation arrives.

Here is a simple implementation in python code.

def auto_gradient(
        f: DifferentiableFunc: Callable[[np.ndarray], np.ndarray],
        args: Sequence[np.ndarray],
        eps: float = 1e-7,
    ) -> Sequence[np.ndarray]:
    """
    Compute the gradient of f with respect to each argument in args.
    """
    grads = []
    for i, arg in enumerate(args):
        grad = np.empty_like(arg) # allocate space for the gradient
        for ix in itertools.product(*[range(x) for x in arg.shape]): # iterate over all indices
            v = arg[ix].copy() # save the original value
            arg[ix] = v + eps
            v1 = float(f(*args))  # must be scalar
            arg[ix] = v - eps
            v2 = float(f(*args))
            arg[ix] = v
            grad[ix] = (v1-v2) / (2*eps) # the slope
        grads.append(grad)
    return grads

Unfortunately, this approach for autogradient is not good enough and requires big resources, so it is more advisable to use a computational graph.

Disclaimer : 100% written by the human brain.

References :

Predict & Generate Output From Fizz Buzz Problem with Machine Learning Approach; Deep FizzBuzzNet.

Rama Reksotinoyo — Wed, 14 Sep 2022 06:29:33 +0000

Wait, because the title is really fancy in my head, i will consider using that title for my thesis if I take master studies.

Maybe you are familiar with FizzBuzzProblem, like al-fatihah in Islam, fizzbuzz problems are also mandatory in the world of job application if you are a software engineer. Usually we will be asked to make a loop i in range(n),

if i is divisible by 3 it will print fizz
if i is divisible by 5 it will print buzz
if i is divisible by 3 and 5 or the number 15 it will prints FizzBuzz.

Implemented in python.

def fizzbuzz(number: int) -> str:
    if number % 15 == 0:
        return "fizzbuzz"
    elif number % 3 == 0:
        return "fizz"
    elif number % 5 == 0:
        return "buzz"
    return str(number)
assert fizzbuzz(1) == "1"
assert fizzbuzz(3) == "fizz"
assert fizzbuzz(5) == "buzz"
assert fizzbuzz(15) == "fizzbuzz"

There are many approaches to solving this problem. From printing manually, this is definitely forbidden for you if you are a software engineer. Until using loop + condition. But there is an approach that I think is really cool, namely the machine learning approach, as I saw on Pydata's youtube video by Joel Gruss
So I decided to build it using another deep learning library, namely Pytorch.

Generating Datasets

First step that i wants to do is, generating the datasets by building a function which accepts two parameters, i and numbers of digits which produces output list of encoded numbers into binary.

def binary_encode(i: int, num_digits: int) -> List[int]:
    return np.array([i >> d & 1 for d in range(num_digits)])

I will tests this function with 10 times.

Next, i will encode of the fizzbuzz output, depends on which output supposed to do.

def fizz_buzz_encode(i: int) -> int:
    if   i % 15 == 0: return 3
    elif i % 5  == 0: return 2
    elif i % 3  == 0: return 1
    else: return 0

def fizz_buzz(i: int, prediction: int) -> str:
    return [str(i), "fizz", "buzz", "fizzbuzz"][prediction]

Spliting datasets into training and testing

X = np.array([binary_encode(i, NUM_DIGITS) for i in range(101, 2 ** NUM_DIGITS)])
y = np.array([fizz_buzz_encode(i) for i in range(101, 2 ** NUM_DIGITS)])

X_train = X[100:]
y_train = y[100:]
X_test= X[:100]
y_test= y[:100]

here's the sizes of splitted datasets

Modelling

This was the most interesting process for me, because i can play with hyperpatameters tunning, here's code :

class FizzBuzz(nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.first = nn.Linear(10, 100)
        self.relu = nn.ReLU()
        self.bactnorm1d = nn.BatchNorm1d(100)
        self.output = nn.Linear(100, 4)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        a = self.first(x)
        relu = self.relu(a)
        batcnorm = self.bactnorm1d(relu)
        return self.output(batcnorm)

I picked cross entropy loss as my criterion and Adam as an optimizers, with 1-e3 for learning rate.

At the end of the training process, this is the loss of my last epoch, i set 1000 epochs by the way.
Validation loss: tensor(0.5942, grad_fn=) Training loss: tensor(0.3454, grad_fn=)

Testing

This was pretty good, the machine learning model working correctly at leasy for 100 testing data.

Full code on my github repo

Reference:
Jeol Gruss's Blog

Youtube

Blog

Attest Is All You Need; Goldbach Conjecture.

Rama Reksotinoyo — Sat, 16 Jul 2022 08:24:46 +0000

First of all, the title was inspired by the original paper of transformer model, Attention Is All You Need

I ever thought that i can getting rich easily just to solve an unsolved math problem, but after i surfing the internet with the keyword “Simplest unsolved math problem” I undo my intention because I am stupid in mathematics.

One of the problems that interests me is to try the "Goldbach Conjecture", it is very simple only integers greater than 2 can be added from 2 prime numbers. It's interesting for me not to answer, but to prove it using my computer science background.

Goldbach conjecture was first mentioned by Christian Goldbach in his letter to Euler in 1742. In his letter, Goldbach reported that even numbers greater than or equal to 4 could be written as the sum of two prime numbers, but he failed to prove the truth of his conjecture.

I proved by using the C + + programming language, this is more or less:

First I declare a function that returns a integer numbers.

int is_prime(int n){
    int num = 1;
    for (int i = 2; i < n/2; i++){
        if((n%i) == 0){
            return num-1;
        }
    }
    return num;
}

The next step is create a procedure for the logic of this mathematical problem.

void solve(){
    int number = 0;
    int primes[100000] = {2};
    int j = 0;
    while(true){
        std::cout<<"Input even number:";
        std::cin>>number;
        if(number > 2 && number % 2==0){

            if(primes[j]<number){
                for (int i = primes[j]+1; i < number; i++){
                    if(is_prime(i) == 1){
                        j++;
                        primes[j] = i;
                    }
                }       
            }
            for (int i = 0; i < j; i++){
                for (int k = 0; k < j; k++){
                    if(primes[i] + primes[k] == number){
                        std::cout<<primes[i]<<" + "<<primes[k]<<" = "<<number<<std::endl;
                        break;
                    }
                }
            }
        }
        else{
            std::cout<<"Input must a even numbers"<<std::endl;
        }
    }
}

The last step is just declared the procedure to main function.

int main(){
    solve();
    return 0;
}

Here is the output, let's say i entered the number 120.

Maybe you passed the incorrect grammar or tenses, i'm sorry, i'm bad in english :)

Reference
Wikipedia

The one who discussed it first

List of unsolved math problem