Mosaic multiple temporal large Geotiff files using GDAL with median value at overlapping pixels

Question

I would like to mosaic multi-temporal large geotiff images (multi-bands) with the code below. However, it returns the last pixel ("last write wins") from the input order. How to get median, or mean values from all inputs rather than from only the lastest input at overlapping areas.

import glob
from osgeo import gdal

def mosaic(input_files, output_file, aoi=None):
    # Create an empty VRT file
    vrt_options = gdal.BuildVRTOptions(resolution='highest')
    vrt_ds = gdal.BuildVRT('/vsimem/mosaic.vrt', input_files, options=vrt_options)

    # Create the mosaic as a GeoTIFF
    gdal.Translate(output_file, vrt_ds, format='GTiff')


    # Cut the output by shapefile
    if aoi:
        warp_options = gdal.WarpOptions(cutlineDSName=aoi, cropToCutline=True)
        gdal.Warp(output_file, vrt_ds, options= warp_options)

    # Clean up the temporary VRT file from memory
    vrt_ds = None
    gdal.Unlink('/vsimem/mosaic.vrt')

I intend to use GDAL and vrt to better process large geotiff files. The corresonding code gives lastest values at overlapping pixels. I expect to get the mean or median values from all inputs in the final output.